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ABSTRACT 

This  report  presents  an  investigation  of  the  effects  that  may  have  an  influence  on  the 
development  of  a  linear  stress  transfer  function  (STF)  relating  the  stress  in  dynamic 
components  to  the  stress  in  static  components.  Effects  such  as  buckling,  non¬ 
uniqueness,  vibration,  and  solution  procedure  are  considered.  Two  procedures  for 
determing  the  STF  are  compared,  one  termed  the  vector  procedure  and  the  other  the 
matrix  procedure.  A  simple  two  dimensional  truss,  which  models  an  idealised 
helicopter  structure,  is  constructed  to  numerically  simulate  the  development  of  a  STF. 
Using  random  inputs  the  resulting  stresses  are  evaluated  exactly.  Noise  is  then  added 
to  both  the  input  loads  and  output  stresses  to  develop  a  noisy  data  set.  Using  this  noisy 
data  set,  STFs  are  developed  using  both  the  vector  and  matrix  techniques.  The  vector 
procedure  is  shown  to  be  sensitive  to  collinearity  in  the  input,  while  the  matrix 
technique  is  found  to  be  more  stable  under  the  same  ill-conditioning. 
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Executive  Summary 

At  present,  helicopter  component  fatigue  estimates  are  necessarily  conservative  to  take 
account  of  the  large  variability  in  component  load  under  the  same  flight  condition.  This 
means  that  most  components  could  still  be  safely  used  past  their  component  retirement 
times.  However,  one  researcher  has  found  that  the  level  of  conservatism  applied  may 
not  always  be  adequate,  and  that  the  retirement  lives  of  a  few  components  may  need  to 
be  reduced  by  as  much  as  75%.  Hence,  some  components  may  remain  in  service  for 
periods  in  excess  of  their  safe  fatigue  lives.  An  accurate  history  of  component  loads 
would  fulfil  two  seemingly  contradictory  objectives,  increasing  safety  while  reducing 
costs. 

However,  due  to  a  harsh  working  environment,  strain  measuring  devices  placed 
directly  on  dynamic  components  have  a  limited  life  span.  In  addition,  the  information 
transmitted  (either  via  slip-rings  or  telemetry)  from  these  in  situ  devices  is  often 
unreliable.  Hence  we  developed  a  linear  stress  transfer  function  (STF)  that  relates  the 
stress  in  a  dynamic  component  to  the  stresses  in  static  components.  An  investigation  of 
STF  development  was  undertaken  using  an  idealisation  of  a  helicopter's  structure. 

Potential  applications  of  this  work  include: 

•  the  accurate  determination  of  component  fatigue  life  expended  (with  component 

retirement  life  implications), 

•  the  parallel  development  of  a  health  monitoring  system  (using  redundant 

information  from  the  STF),  and 

•  a  real-time  system  warning  pilots  of  excessive  aircraft  loads. 
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1.  Introduction 


Helicopter  fatigue  damage  estimates  are  necessarily  over-conservative  at  present,  and 
do  not  account  for  variations  in  usage  on  an  aircraft-to-aircraft  basis.  This  can  result  in 
maintenance  intervals  that  are  shorter  than  necessary  and  the  premature  replacement 
of  some  components  [1].  As  Gunsallus  and  Robeson  [2]  correctly  point  out,  '[t]he  user  is 
constantly  pushing  the  envelope  of  safe  operation  in  efforts  to  extract  the  most  safe  usage  from  a 
given  airframe'.  This  explains  the  emphasis  on  improving  the  estimation  of  helicopter 
fatigue  damage  [3-17]. 

Helicopters,  unlike  fixed-wing  aircraft,  usually  have  single  (non-redundant)  load  paths 
and  are  subjected  to  both  high  cycle  and  low  cycle  fatigue  loads  [18].  The  fatigue  life  of 
a  helicopter  component  is  clearly  dependent  on  the  loading  history  experienced  by  that 
component,  but  the  manufacturer's  fatigue  life  estimation  is  based  on  a  predicted 
usage  spectrum  [1, 19-25],  However,  the  usage  spectrum  experienced  by  helicopters  is 
normally  different  to  that  envisaged  by  the  operator  and  the  manufacturer  before  the 
aircraft  enters  service.  This  varied  usage  spectrum,  together  with  the  fact  that  the 
severity  of  individual  manoeuvres  varies  widely  both  between  pilots  and  between 
individual  helicopters  [26,  27],  imply  that  the  fatigue  lives  of  helicopter  components  are 
necessarily  overly  conservative  in  nature.  An  accurate  estimation  of  loads  in  critical 
components  would  have  two  seemingly  contradictory  benefits  —  improved  safety  and 
reduced  component  replacement  costs.  There  are  possibly  two  added  bonuses  of 
installing  an  effective  on-line  load  detection  system.  The  first  is  a  real-time  warning 
system  that  could  alert  the  pilot  that  a  high  fatigue  manoeuvre  is  being  entered,  thus 
allowing  the  incidence  of  damaging  loads  to  be  reduced  [28],  The  second  is  a  fault 
warning  capability  (either  real-time  or  off-line),  which  could  be  implemented  using 
redundant  information  from  the  developed  system. 

Currently,  after  every  flight,  the  helicopter  aircrew  of  the  Australian  Defence  Force  has 
to  complete  an  EE360  Fatigue  Monitoring  Form,  which  documents  rudimentary  and 
easily  quantifiable  data  on  manoeuvres  etc.  For  example,  the  information  recorded  on 
EE360  forms  includes  the  number  of  landings,  type  of  landing  surface,  percentage  of 
time  at  high  and  low  altitudes,  and  stores  configuration.  As  expected,  only  simple 
fatigue  calculations  can  be  performed  using  EE360  data. 

Barndt  and  Moon  [19,  p.  1375]  state  that  installing  a  Structural  Data  Recording  Set 
(SDRS)  system  ‘translates  into  millions  of  dollars  saved  even  after  the  installed  system  costs 
are  considered'  and  that  a  SDRS  system  would  help  ' identify  aircraft  exposed  to  severe 
operating  conditions'.  Fatigue  life  estimates  based  on  regime  recognition  provide  an 
improvement  when  compared  to  those  based  solely  on  usage  hours,  but  they  still  do 
not  fully  account  for  the  actual  variation  of  loads  during  a  given  manoeuvre  [29]. 

The  next  level  of  complexity  involves  measuring  stresses  directly  on  critical  structural 
components,  which  are  mostly  from  the  rotor  system  for  helicopters.  (For  example,  of 


1 


DSTO-KR-0171 


the  21  most  fatigue-critical  components  on  the  AH-64A  Apache,  Lombardo  [25,  p.  46] 
states  that  20  are  from  the  rotor  system.)  However,  slip  rings  or  rotating  telemetry 
required  for  such  direct  measurements  of  stress  lack  the  reliability  and  maintainability 
needed  for  everyday  fleet  usage  [29].  This  constraint  forces  us  to  seek  an  alternative 
approach,  such  as  developing  a  stress  transfer  function  that  relates  stresses  experienced 
by  rotor  system  components  to  stresses  measured  in  fixed  components. 

A  companion  report  [30]  surveys  published  research  on  helicopter  fatigue  estimation 
techniques.  However,  there  appears  to  be  a  lack  of  elementary  research  into  the 
underlying  structure  of  this  estimation  problem.  For  example,  several  papers  tackle  the 
problem  of  non-orthogonal  inputs  using  trial-and-error.  Although  somewhat  effective, 
these  procedures  may  not  lend  themselves  to  generalisations,  especially  when  different 
helicopter  types  are  considered.  In  this  report,  we  investigate  some  of  the  more 
fundamental  questions  concerning  the  development  of  a  stress  transfer  function. 
Although  the  development  is  not  exhaustive,  the  findings  do  highlight  some 
interesting  properties  about  the  underlying  structure  of  the  problem. 

Throughout  this  report  we  assume  the  outputs  from  strain  gauges  have  been  converted 
to  stresses.  As  such  we  use  strain  gauge  output  and  stress  inter-changeably.  We 
investigate  aspects  that  influence  the  development  of  stress  transfer  functions  between 
strain  gauge  outputs  in  static  components.  (We  envisage  that  in  the  future  this  work 
will  be  expanded  to  include  dynamic  effects.)  The  effects  investigated  in  developing 
transfer  functions  include,  non-linearities  (in  particular  those  due  to  buckling),  non¬ 
uniqueness  (for  example,  rotor  blade  stress  distribution),  and  solution  procedure  (how 
loading  information  affects  results). 

A  simplified  statically  determinate  helicopter  structure  is  developed  to  examine 
alternative  techniques  of  setting  up  the  transfer  functions.  The  difference  between  these 
techniques  is  whether  the  externally  applied  loads  are  known  or  unknown,  and  due  to 
the  nature  of  the  solutions  we  term  them  matrix  and  vector  techniques  respectively. 
Either  technique  may  be  used  for  rotor  component  stress  estimation,  however,  the 
matrix  technique  requires  knowledge  of  external  loading  for  the  development  of  the 
transfer  functions.  Once  the  transfer  functions  are  developed,  neither  technique 
requires  external  loading  information.  When  the  stress  measuring  systems  contain  no 
noise  the  two  techniques  give  rise  to  identical  solutions.  However,  initial  results 
suggest  that  the  matrix  solution  technique  is  more  robust  when  the  underlying  system 
of  equations  is  ill-conditioned.  In  fact,  when  the  condition  number  (a  measure  of 
sensitivity)  of  the  underlying  system  of  equations  is  large,  the  vector  solution  technique 
may  be  unusable  even  for  small  amounts  of  noise  in  the  stress  measurement  system. 

The  appendices  present  the  investigation  of  a  number  of  problems  that  arose  along  the 
way.  Appendix  A  shows  that  indeterminate  systems  are  still  linear.  In  Appendix  B  we 
analytically  solve  for  the  stress  at  eight  locations  in  a  simplified  helicopter  truss.  The 
truss  is  two-dimensional  and  statically  determinate.  Appendix  C  examines  the  effect  of 
a  static  assumption  on  a  dynamically-loaded  simply-supported  beam.  Conditions 
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under  which  the  static  assumption  remains  valid  are  stated.  Appendix  D  gives  a  brief 
overview  of  the  least  squares  solution  from  a  singular  value  decomposition 
perspective.  Finally,  Appendix  E  proves  that  the  error  surface  (as  defined)  between  a 
matrix  and  its  approximation  has  only  two  maxima  and  two  minima. 

Applications  of  this  work  include: 

•  determining  the  fraction  of  fatigue  life  expended  by  components, 

•  a  health  monitoring  system  (for  example,  a  system  which  detects  cracks),  and 

•  a  real-time  warning  system  of  excessive  loads. 
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2.  Truss  Model 


By  developing  a  simple  two-dimensional  analytic  model  in  this  section  we  can  readily 
gain  information  about  what  problems  may  arise  in  more  complex  models  (including 
the  actual  helicopter). 


2.1  Simple  Truss  Model 

We  first  consider  the  highly  simplified  truss  model  shown  in  Figure  2.1  for  a  generic 
helicopter.  The  striped  sections  above  the  main  and  tail  rotors  represent  loading,  wad 
and  niAC  are  point  masses,  G,  for  i=0,l,2,...,7  denote  strain  gauges,  Px,  P}/,  and  Pm  are  the 
main  rotor  loads,  and  TX/  Ty,  and  Tm  are  the  tail  rotor  loads.  All  members  are  pinned 
together  except  for  the  vertical  member  connected  to  the  main  rotor,  which  passes 
through  a  slide  bearing.  The  point  masses  allow  forces  due  to  accelerations  (gravity  for 
example)  to  be  simulated. 
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As  can  be  seen  from  Figure  2.1  all  main  rotor  system  loads  are  transferred  into  the 
fuselage  via  the  main  rotor  shaft  (we  will  see  in  §  2.3  that  this  causes  uniqueness 
problems  when  determining  blade  stresses  from  fuselage  stresses).  The  tail  rotor  may 
be  idealised  in  the  same  way.  For  simplicity  we  configured  the  problem  to  be  a  static 
one  by  pinning  the  bottom  left  joint  to  the  ground  and  adding  a  roller  to  the  bottom 
right  joint  (located  to  the  right  of  gauge  G2).  Appendix  A  shows  two  sets  of  examples 
demonstrating  why  indeterminate  systems,  although  more  complicated,  may  be 
treated  using  the  same  linear  principles.  The  solution  for  the  statically  determinate 
truss  shown  in  Figure  2.1  is  given  in  Appendix  B.  Vibration  effects  in  relation  to  the 
static  assumption  are  discussed  in  Appendix  C,  where  we  show  that  under  suitable 
conditions  the  static  assumption  yields  a  good  approximation. 


2.2  Introduction  of  Non-linearity 

One  way  that  the  simple  truss  shown  in  Figure  2.1  could  become  non-linear  is  if  one  of 
the  components  were  to  buckle.  (On  an  actual  helicopter  the  skin  is  an  example  of  a 
component  subject  to  buckling.)  Consider  the  rectangular  shear  panel,  shown  in  both 
unbuckled  and  buckled  states,  in  Figure  2.2  (figure  after  Brunh  [31,  p.  C11.8]  and  in 
Timoshenko  and  Gere  [32,  p.  421]).  We  will  consider  what  effects  skin  buckling  has  on 
the  supporting  stiffeners. 


Figure  2.2.  Buckled  shear  panel  with  ellipses  depicting  wrinkling.  The  shear  panel 
has  a  bi-linear  stress  response  to  the  vertical  load,  linear  both  before  and  after 
buckling. 

The  stresses  on  the  left  side  of  the  panels,  fc  and  ft ,  represent  compressive  and  tensile 
stresses  respectively.  The  shear  panel  on  the  left,  under  the  critical  load  Vcr ,  is  in  a  state 
just  prior  to  buckling.  The  shear  panel  on  the  right,  under  a  load  Vcr  +  Vt ,  has  buckled 
showing  the  characteristic  wrinkling  (depicted  as  ellipses).  After  buckling,  the  shear 
panel  on  the  right  may  be  thought  of  as  being  composed  of  two  shear  panels:  one  panel 
just  before  buckling  (left  panel)  that  can  take  compression,  and  another  panel  that  can 
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act  only  as  a  tension  field  (middle  panel).  As  can  readily  be  seen  any  stress  reading 
taken  as  a  function  of  vertical  load  would  result  in  a  continuous  bi-linear  function  with 
different  slopes  before  and  after  buckling.  Before  buckling  the  panel  would  linearly 
take  both  compression  and  tension.  After  a  critical  load  Vcr  however,  the  panel  would 
continue  to  take  only  tensile  load  .  More  precisely,  the  stress  in  the  connecting  stiffener 
would  be  of  the  form 


[  cxV,  F<VCf. 

a  [c2  +  c3V ,  V>Vcr 

where  V  is  the  vertical  load  and  clr  c2,  andc3  are  constants  which  depend  on  the 
shear  panel  and  stiffener  properties  (such  as  stiffness,  area,  and  length)  and  on  the 
angle  which  the  wrinkles  make  with  the  horizontal.  (See  Bruhn  [31,  §01.5]  for  further 
details.) 

This  type  of  non-linearity  applies  to  local  buckling  as  well,  for  example  the  flanges  of  a 
beam.  Thus  the  stress  within  the  beam  will  exhibit  a  bi-linear  response  to  the  applied 
load  if  the  flange  buckles  during  loading.  These  non-linear  effects  need  to  be  taken  into 
account  in  an  implementation  of  this  stress  correlation  method,  and  thus  gauges 
should  be  placed  away  from  buckling  structures  to  achieve  a  linear  relationship 
between  measured  stresses. 


2.3  Problems  with  Uniqueness 

Consider  the  truss  shown  in  Figure  2.3  under  two  loads  Fj  and  F?  at  a  distance  xi  and  X2 
from  the  left  end  of  the  beam  respectively  (the  shaded  beam  has  length  L). 


Figure  2.3.  A  non-unique  loading  configuration. 
Solving  for  the  reaction  forces  Ra,  Rb,  and  Rc  on  the  shaded  beam  gives 
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R 

R 


--R 


B  cos  6a 
cos  dA 


sin(0A 


+ 


(P-a), 


Rc  -  a, 


and 


where  a-F1(x1/L)  +  F2(x2/L )  and  P  =F1+F2.  Thus  the  reaction  forces  are  only  a 
function  of  a  and  p.  So  that  for  any  combination  of  constants  a  and  j 8  the  reactions  are 
also  constant,  despite  the  fact  that  both  the  loads  and  where  they  are  applied  may  be 
varying.  For  example,  if  L=1  then  the  combination  Fj=3,  F2=7,  xj=0.3,  and  X2=0.5  gives 
the  same  stress  distribution  as  the  combination  Fi= 18,  F2=-8,  Xi= 0.6,  and  X2=0.8 
everywhere  on  the  truss  except  in  the  shaded  beam.  There  is  no  unique  way  to 
determine  the  stress  in  the  shaded  beam  from  the  stress  in  the  remainder  of  the  truss 
using  static  information  alone. 

The  same  non-uniqueness  problem  arises  when  the  blade  loads  are  transferred  through 
the  main  rotor  mast  onto  the  underlying  structure.  Figure  2.4  shows  two  different 
loading  distributions  on  a  single  blade  that  produce  the  same  vertical  shear  force  and 
bending  moment  at  the  hub.  (If  the  blade  has  unit  length  then  the  two  bending  moment 
distributions  shown  in  Figure  2.4  are  given  approximately  by  5.486  -  8.869x  +  2.5x2  - 
1.647x3+  3.839x4-  1.310x5  and  5.486-  8.869x  -  34.05x2  +  95.51x3-  80.82x4+  22.74x5. 
Differentiating  the  bending  moment  once  gives  the  vertical  shear  force  and 
differentiating  twice  gives  the  loading.)  As  can  be  seen  from  Figure  2.4,  from  statics 
alone  there  is  no  way  to  determine  what  the  blade  stress  distribution  is  from  the 
structure  below  the  main  rotor  hub,  since  the  shear  and  bending  moment  are  identical 
at  the  hub. 


Load  Shear  Moment 


Figure  2.4.  Tire  hub  sees  the  same  shear  and  bending  moment  for  two  different 
loading  distributions. 


It  is  easy  to  prove  that  there  is  an  infinite  number  of  loading  distributions  that  will  give 
rise  to  any  single  combination  of  shear  and  bending  moment.  Let  the  bending  moment 
distribution  be  given  by 


M(x)  =  a  +  fx  +  /(x) , 
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where  a  and  ft  are  constants,  and/(x)  is  any  function  that  satisfies  the  two 
conditions  /(0)=  /'(0)  =  0 .  Since  shear,  V,  and  the  loading  distribution,  F,  are  simply 
given  by  the  first  and  second  derivatives  of  bending  moment  respectively,  we  have  that 

V(x)=  p  +  f'(x), 

and 

f(*)  =/'(*)• 

Thus  the  bending  moment  and  shear  at  the  origin  will  always  be  M(0)  =  a  and 
V(0)=  P  respectively,  regardless  of  the  function  f(x)  chosen  (provided  it  and  its  first 
derivative  are  zero  at  the  origin). 

Although  in  theory  there  is  an  infinite  number  of  functions  that  satisfy  the  same 
bending  and  shear  values  at  the  root,  in  practice  a  vast  majority  of  these  loading 
distributions  would  never  be  realised.  However,  there  are  still  a  sufficiently  large 
number  of  possible  load  combinations  to  yield  a  difficult  non-uniqueness  problem. 


2.4  Stress  Estimation  with  Known  External  Load 

In  the  remainder  of  this  report  we  assume  all  structures  are  static,  that  is,  we  ignore  all 
dynamic  effects.  Define  a  load  as  any  input  having  an  influence  on  stress  so  that,  for 
example,  shearing  forces  (such  as  Px),  bending  moments  (such  as  Tm),  and  accelerations 
(such  as  gx)  are  defined  as  loads.  If  we  know  both  the  loading  and  resulting  stress  on 
our  simple  helicopter  truss  we  can  write  down  the  relation  for  the  stress  at  the  location 
of  the  kth.  strain  gauge  as  a  function  of  load 

G k  (0  =  aikJ\  (0  ■*"  aik  /a(0  ^  h  ajkf  j  (0  '  (2-1) 

where  the  ajk  are  constants  (the  subscripts  j  and  k  refer  to  the  load  and  stress 
respectively),  /  is  the  number  of  input  loads,  and  /.  is  the  ;th  load.  (We  have 
implicitly  assumed  that  stress  is  a  linear  function  of  input  loads.)  Introduce  a  bracketed 
superscript  notation  for  time  so  that  <7^  =crk(ti) ,  and  apply  a  similar  notation  for  the 
loading  parameters.  We  can  then  write  a  linear  system  of  equations  that  relate  stress  to 
loading  as 
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where  I  is  the  number  of  observations  (of  load  and  stress  taken  at  different  times)  and 
K  is  the  number  of  strain  gauges. 

For  ease  of  description  let's  consider  a  single  strain  gauge,  that  is  K  =  1 .  In  this  case  the 
system  of  equations  is  over-determined  if  I  >  J  (there  are  more  observations  I  than 
unknowns  J )  and  under-determined  if  I  <  J  (there  are  more  unknowns  than 
observations).  In  practice  the  under-determined  system  would  not  arise  if  an  arbitrary 
number  of  observations  were  possible,  and  hence  will  not  be  discussed  below. 

If  the  system  has  full  rank  (that  is,  the  force  matrix  has  linearly  independent  rows)  and 
there  are  as  many  observation  as  loads  ( I  =  J )  then  we  can  easily  find  the  coefficient 
matrix.  The  product  of  the  inverted  force  matrix  by  the  stress  matrix  gives  the 
coefficient  matrix 

[AHfm 

where  [A]  is  the  force  coefficient  matrix,  [F]  is  the  force  matrix,  and  [S]  is  the  stress 
matrix.  However,  if  we  were  to  have  the  luxury  of  being  able  to  choose  the  force  matrix 
to  be  the  identity  matrix,  then  the  coefficient  matrix  would  be  identical  to  the  stress 
matrix. 

If  the  system  has  full  rank  and  is  over-determined  (/>/),  then  there  is  no  unique 
solution.  A  unique  solution  is  obtained  by  minimising  the  norm  of  the  residual,  where 
the  residual  is  defined  as  [FlA]-[Z] .  Minimising  over  the  2-norm  gives  rise  to  the  least 
squares  (LS)  problem.  Golub  and  van  Loan  [33]  suggest  either  the  method  of  normal 
equations  or  QR  factorisation  for  the  solution  of  the  LS  problems.  The  method  of 
normal  equations  involves  computing  the  Cholesky  factorisation  of  the  product  of  the 
force  matrix  transpose  and  the  force  matrix.  The  accuracy  of  this  method's  solution 
depends  on  the  square  of  the  condition  number  K ,  which  for  the  2-norm  is  defined  as 
the  ratio  of  the  largest  to  smallest  singular  values  of  the  force  matrix  (see  Appendix  D.l 
for  a  discussion  of  singular  values).  The  QR  factorisation's  LS  solution  involves 
obtaining  the  Householder  QR  factorisation  of  the  loads  matrix  and  using  the  result  in 
the  solution  of  an  equivalent,  but  easier  to  solve,  LS  problem.  The  QR  factorisation 
method  runs  into  trouble  whenever  k~  1/u  (where  u  is  the  computer's  unit  round 
off  error,  that  is  1  +  u  =  1 ),  while  the  method  of  normal  equations  is  less  robust  running 
into  trouble  whenever  k~  1/  Vu  .  The  two  methods  produce  comparable  inaccuracies 
when  applied  to  large  residual,  ill-conditioned  problems.  Golub  and  van  Loan  also 
demonstrate  that  the  condition  number  grows  if  a  column  is  added  to  the  force  matrix 
(that  is,  an  extra  load  is  added  to  the  system). 

If  the  system  of  equations  is  rank  deficient  (that  is,  at  least  one  of  the  rows  from  the 
force  matrix  is  a  linear  combination  of  the  remaining  rows),  then  the  LS  problem  has  an 
infinite  number  of  solutions.  Singular  value  decomposition  (SVD)  is  particularly  useful 
for  rank  deficient  problems  since  it  provides  a  neat  expression  for  the  LS  solution  and 
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the  norm  of  the  minimum  residual  [33,  p.  242].  (Appendix  D.l  contains  details  of  the 
least  squares  SVD  solution.) 

From  Equation  (2.1)  the  stress  at  the  zeroth  gauge,  in  terms  of  loads,  is  given  by 
e0={A0Y{f}+{B0Y\gl  where  {A0}  and  {B0}  are  the  load  and  acceleration  coefficient 
vectors  for  the  zeroth  stress,  and  {/}  and  {g}  are  the  force  and  acceleration  vectors.  (Note 
that  the  load  vector  has  been  partitioned  into  a  force  and  acceleration  vector,  since  we 
will  in  general  know  the  acceleration  vector.  In  contrast  we  will  not  know  the  force 
vector,  at  least  not  after  the  force  coefficient  matrix  correlation  is  set  up.)  Similar 
expressions  apply  to  the  remaining  stresses  so  that 

M=  [A]{/}+  [B]{g}, 

where  curly  and  square  brackets  denote  vector  and  matrices  respectively, 
{c}={(T1  <r2  <7 3  C74  c5  cr6Jr  is  the  stress  vector,  [A]e916><6  and  [B]e9l6><2  are  the 

force  and  acceleration  coefficient  matrices  respectively,  and 
{/}  =  (/i  /2  h  fi  /5}T  and  }g}  =  {g,  g2V  are  the  force  and  acceleration  vectors 

respectively.  Initially  we  will  use  only  six  of  the  seven  stress  equations  so  that  the  force 
coefficient  matrix  A  is  square.  Also  note  that  the  stress  vector  cr  has  components  Oi, 
02,...,  06,  which  denote  any  six  stress  equations,  not  necessarily  the  first  six.  Assuming 
the  force  coefficient  matrix  is  invertible,  solving  for  the  force  vector  in  terms  of  the 
stress  vector  gives 


{fHAHM-Mk]). 

And  hence  the  stress  at  the  zeroth  gauge  in  terms  of  the  stress  vector  is 

<6,  ={A,}r[Ar({^}-[BW+{B,7fe}- 

The  accelerations  need  not  have  been  partitioned  as  they  have  been  above;  we  could 
have  left  them  in  the  loading  matrix.  However,  because  we  can  normally  easily  obtain 
the  accelerations,  we  might  as  well  subtract  them  from  the  system  we  are  trying  to 
solve  (as  was  done  in  the  above  expression).  The  above  formulation  should  result  in  a 
more  stable  system  than  would  arise  if  the  accelerations  had  been  left  in  the  loading 
matrix  as  unknowns. 

To  make  the  following  procedures  clearer  let's  ignore  the  acceleration.  We  may  then  re¬ 
write  the  zeroth  stress  as 


e0={A0}T[AF{o}. 


(2.2) 
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2.5  Stress  Estimation  with  Unknown  External  Load 

If  all  we  know  are  the  stresses  at  a  certain  number  of  locations  then  it  is  possible  to 
determine  the  stresses  at  still  other  locations.  More  specifically,  if  we  know  the  stresses 
at  K  strain  gauge  locations  then  the  stress  at  a  different  location,  say  the  zeroth  location, 
is  given  by 


CTq  —  OC^CJ]  +  CH^O 2  +  *  *  *  +  CCg(J fr 

=  {af{al 


(2.3) 


where  the  ak  are  constants,  and  {a}  and  {a}  are  the  vectors  containing  the  constant  ak 
and  stresses  ak  respectively.  Such  a  relation  exists  if  there  are  as  many  strain  gauges  as 
there  are  input  loads  and  the  stresses  on  the  left  are  of  full  rank.  That  is,  no  stress  on 
the  right  hand  side  (RHS)  is  a  linear  combination  of  other  stresses  on  the  RHS,  or  in 

mathematical  notation  o, ■  *  ]Tf=i Pkok  for  any  i  =  1,2,..., K  ,  where  the  (3k  are  constants. 


We  now  see  why  we  termed  the  solution  technique  involving  the  load  measurement 
the  matrix  technique,  since  Equation  (2.2)  involves  a  matrix.  On  the  other  hand, 
Equation  (2.3),  developed  without  external  load  knowledge,  and  only  involving  a 
vector,  is  termed  the  vector  technique.  The  other  difference  to  note  is  that  we  can  think 
of  the  matrix  technique  as  an  indirect  method,  but  the  vector  technique  as  a  direct 
method.  After  all,  the  matrix  and  vector  techniques  yield  the  same  results  under  ideal 
conditions  (no  noise),  and  so  developing  the  force  coefficient  matrix  seems  an 
unnecessary  step  in  the  matrix  technique  for  producing  the  same  output  as  the  vector 
technique. 
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3.  Numerical  Example  of  a  2-D  Helicopter  Truss 

In  this  section  we  numerically  simulate  stresses  (at  strain  gauge  locations)  on  a  2-D 
helicopter  truss  (using  the  information  derived  in  Appendix  B)  to  determine  the 
stability  properties  of  developing  a  stress  transfer  function.  All  numerical  values  given 
in  this  section  have  been  non-dimensionalised.  The  pin-joint  locations  of  a  simple  truss 
(see  Figure  B.l)  are  shown  in  Table  3.1,  while  Table  3.2  shows  the  location  of  the  two 
point  masses  and  eight  strain  gauge  locations. 


A 

B 

C 

D 

E 

F 

G 

H 

X 

0 

3 

2 

2 

2 

7.5 

3.6 

3.3 

y 

0 

0 

1 

0 

1.5 

0.8 

0.6 

0.9 

Table  3.1  Location  of  helicopter  truss  pin-joints. 


P 

niAc 

mA  e 

Go 

Gi 

g2 

g3 

g4 

g5 

g6 

g7 

X 

0.5000 

0.5000 

2.0000 

1.0000 

2.5000 

2.6667 

1.5000 

2.4333 

3.0750 

3.4500 

lx 

0.2500 

0.0000 

0,5000 

0.5000 

0.0000 

0.3333 

0.0000 

0.9667 

0.2250 

0.4500 

Table  3.2  Location  of  helicopter  truss  masses  and  strain  gauges. 


The  two  masses  were  mAc= 9  and  mAE=7,  and  all  beam  members  were  circular  tubes 
(outer  radius  0.015  and  thickness  0.005)  with  a  cross-sectional  area  of  392.7x10-6  and  a 
second  moment  of  area  of  31.91xl0-9.  Using  these  physical  properties  the  stress  at  the 
eight  strain  gauge  locations  (shown  below  with  five  significant  figures)  are 

c 70  =  -235.06  x  103  Pm  + 117.53  x  103  Px  +  2.5465  x  103  Py , 

Gt  =  -1.8980  xlO3  Pm  +  2.8471xl03  Px  -  529.74  xlO3  gx  + 1.0595  xlO6  gy 
- 1.8980 xlO3  Tm  +1.5184 xlO3  Tx  -8.5412 xl03Ty, 

<y2  =  -848.83  Pm  + 1.2732 xlO3  Px  + 156.71  x  103  Py  +1.9099 xlO3  gx 

+  270.42 xlO3  gy  +  1.6977xl03  Tm  +1.1884xl03  Tx  +7.6394 xl03Ty, 

<73=1.2004x103  Pm  -1.8006 xlO3  Px  - 2.7009 xlO3  gx  +  5.4019 xlO3  gy 
+ 1.5005 xlO3  Tm  -1.2004 xlO3  Tx  +6.7524x10 3Ty, 


(3.1) 

(3.2) 

(3.3) 

(3.4) 
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aA  =  1.6977  xlO3  Pm  +  235.06  xlO3  Py  + 1.9099  xlO3  gx  +  818.89  xlO3  gy 
+  1.6977 xlO3  Tm  + 1.1884 xlO3  Tx  +7.6394 xl03Tv, 

(j5  =-2.7668 xlO3  Tm  +2.2135 xlO3 Tx  -12.541  xlO3 Ty/ 

ct6  =  -4.6974  x  103  Tm  -  268.427;  - 17.112  x  103  Ty , 

cr7  =6.0021  xlO3  Tm  +  600.217;  +  25.209 x  103 Ty . 


(3.5) 

(3.6) 

(3.7) 


(3.3) 


Using  either  the  external  loads  or  the  stresses  (from  gauges  1  to  7)  to  set  up  a  coefficient 
matrix  or  vector  respectively,  we  will  develop  a  stress  transfer  function  between  the 
zeroth  gauge  and  six  other  gauges. 

3.1  Noiseless  Stress  Estimations 

In  this  section  we  investigate  what  problems  arise  when  we  try  to  estimate  stress,  using 
measurements  with  no  noise  (either  in  the  stress  measurements  or  the  load 
measurements).  Then  we  investigate  what  effect  measurement  noise  has  on  our  stress 
estimation  models. 

Assuming  no  measurement  noise  we  can  readily  derive  the  stress  equations  given  by 
Equations  (3.1)-(3.8).  A  small  amount  of  work  quickly  shows  that  the  stress  equations 
02,  03,  and  05  are  linearly  dependent;  in  fact,  the  following  relation  holds 

<t,  =-1.5811<t3  -0.17150(7 5 . 

Thus  if  we  eliminate  any  one  of  the  stress  equations  1,  3,  or  5  the  resulting  force 
coefficient  matrix  will  be  invertible,  because  the  system  will  then  have  full  rank.  What 
would  happen  if  we  eliminated  any  other  row  (that  is  not  1,  3,  or  5),  say  for  example  the 
seventh  stress  equation  o>?  Then  the  resulting  force  coefficient  matrix  A  would  be  a 
non-unique  linear  transformation,  with  the  null  space  vector  (explained  below)  given 
by 


{n}  =  { -0.12308  0.00400014  -  0.18462  -  0.24616  -  0.24616  0.91080  }r . 

By  definition  the  matrix  [A]  maps  the  null  space  vector  {n}  onto  the  zero  vector  {0},  that 
is  [A]{n}={0).  Thus  if  [A]{x}  =  {6} ,  then  for  any  scalar  value  a  the  vector  {x}  +  a{n}  has 
the  same  solution  under  the  linear  transformation  A.  In  other  words,  if  { y }  =  {x}  +  a{n} 
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then  the  product  [A]{y}  also  equals  {b},  despite  the  fact  that  (x)  and  {y}  are  different 
vectors. 

If  there  are  more  stress  equations  than  unknowns,  then  singular  value  decomposition 
(SVD),  see  Golub  and  van  Loan  [33,  pp.  246-248],  is  one  way  to  determine  which  rows 
to  eliminate  The  number  of  singular  values  above  some  small  lower  limit  (see 
Appendix  D.l  for  further  details  on  SVD)  gives  the  rank  of  the  matrix.  As  Golub  and 
van  Loan  state:  'Near  rank  deficiency  in  A  cannot  escape  detection  when  the  SVD  of  A 
is  computed'.  The  2-norm  condition  number  k  of  the  matrix  A,  given  by  the  ratio  of 
largest  to  smallest  singular  values  (see  Appendix  D.l),  is  a  convenient  by-product  of 
using  SVD  to  determine  the  rank  of  A. 

To  improve  stability  properties,  we  choose  the  resulting  coefficient  matrix  that  yields 
the  smallest  condition  number.  Sequentially  eliminating  either  the  first,  third,  or  fifth 
stress  equations  (oi,  <J3,  and  a5)  from  the  force  coefficient  matrix  yields  condition 
numbers  x=569.97,  k= 562.97,  and  k=1805.8  respectively.  So  to  minimise  the  condition 
number  we  eliminate  the  third  stress  equation  03,  which  results  in  the  force  coefficient 
matrix 


2847.1 

0 

-1898.0 

1518.4 

-8541.2 

-1898.0' 

1273.2  156700. 

-848.83 

1188.4 

7639.4 

1697.7 

0 

235060. 

1697.7 

1188.4 

7639.4 

1697.7 

0 

0 

0 

2213.47 

-12450.8 

-2766.84 

0 

0 

0 

-268.42 

-17112. 

-4697.4 

0 

0 

0 

600.21 

25209. 

6002.1 

The  columns  represent  the  input  loads  Px,  Py,  Pm,  TXr  Ty,  and  Tm  from  left  to  right 
respectively.  The  above  matrix  has  linearly  independent  rows  and  a  2-norm  condition 

number  k=562.97. 


From  the  zeroth  stress  equation  Co,  given  by  Equation  (3.1),  we  know  that  the  force 
coefficient  vector  is  given  by  A0T  =  103 {117.53  2.5465  -235.06  0  0  0},  and  using 


Equation  (2.2)  we  have  that 


117.53  ' 

T 

‘  614.67 

-589.05 

2.5465 

-2.8538 

6.3814 

-235.06 

395.15 

-883.57 

0.0000 

>• 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

392.70 

-346.36 

62.091 

138.84 ' 

tii 

0.0000 

-0.48943 

-2.0180 

-4.5123 

<72 

589.05 

-158.12 

93.137 

208.26 

ct4 

0.0000 

391.54 

124.18 

277.68 

<?5 

0.0000 

-30.119 

372.55 

277.68 

<?6 

0.0000 

87.344 

-1577.1 

-1027.4 

?1. 
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(Tq  =  —20.6479<t |  +138.478a2  -92.3077ct4  -3.54108<t5  -14.6003ct6  -32.6472ct7  (3.9) 

As  can  be  seen  from  the  above  expression  the  third  stress  equation  03  was  not  used  to 
solve  the  above  system,  thus  yielding  a  system  with  full  rank  and  the  smallest  possible 
condition  number.  Remember  that  the  above  solution  was  produced  with  no  noise  in 
the  system.  We  now  investigate  the  effects  of  noise  in  the  both  load  and  stress 
measurements. 


3.2  Noisy  Stress  Estimations  with  Known  External  Load 
(matrix  technique) 

In  this  report  we  refer  to  randomly  generated  numbers.  All  random  numbers  were 
generated  with  a  pseudo-random  number  generator  defined  within  the  Mathematica 
package  [34],  For  consistency,  all  comparisons  involving  random  numbers  were  carried 
out  with  the  same  set  (or  same  initial  set  if  subsequent  comparisons  involved  larger 
sets)  of  random  points.  The  same  set  of  random  points  was  easily  re-generated  by 
seeding  the  pseudo-random  number  generator  with  the  same  starting  integer. 

Furthermore,  assume  we  have  a  priori  knowledge  of  the  exact  force  coefficient  matrix, 
and  hence  know  that  eliminating  the  third  stress  equation  ct3  leads  to  the  smallest 
condition  number.  In  reality  no  such  a  priori  knowledge  would  be  available,  and  thus 
an  appropriate  way  to  determine  which  combination  of  strain  gauges  minimises  the 
condition  number  for  noisy  data  needs  further  investigation.  However,  we  later  show 
that  for  this  particular  truss,  eliminating  the  third  stress  equation  does  indeed  lead  to 
'optimum'  results. 

Input  loads  were  randomly  generated  uniformly  in  the  range  [-10, +10]  for  each  of  the 
six  input  loads  Px,  Py,  Pm,  Tx,  Tv,  and  T„„  resulting  in  1000  different  load  sets.  Each  set 
contained  six  input  loads,  for  simplicity  the  accelerations  gx  and  g}/  were  set  to  zero. 
These  random  input  loads  (together  with  statics)  were  then  used  to  generate  1000  stress 
sets  (each  set  containing  seven  stresses  from  gauges  0  to  7  excluding  gauge  3).  In 
summary,  we  randomly  selected  the  input  loads,  but  the  resulting  truss  stresses  were 
determined  using  statics.  We  interchangeably  refer  to  a  load  set  or  a  stress  set  as  a  load 
point  or  stress  point  respectively.  White  noise  was  then  added  to  both  the  resulting 
loads  and  stresses,  at  5%  and  10%  of  the  recorded  value.  The  stress  equations  were 
again  derived,  this  time  using  a  least  squares  (LS)  fit  of  these  noisy  data  sets.  (The  in¬ 
built  Mathematica  function  Fit  []  ,  based  on  the  pseudo-inverse  of  the  input  data, 
performs  the  LS  fit.)  The  force  coefficient  matrices  for  5%  and  10%  noise  are  given 
respectively  by 
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■^■5%  ~ 


and 


Ao% - 


2847.9  -20.972 

1209.9  156880. 

-93.037  235310. 
5.3447  -31.073 

-1.1641  -38.491 
0.10014  55.116 


2844.1  -41.703 

1160.9  156790. 

-160.52  235170. 
11.159  -62.219 

-1.2640  -77.541 
-1.4135  111.01 


-1883.3  1515.5 

-947.81  967.17 

1564.1  850.91 

27.757  2204.8 

32.788  -  286.31 

-48.468  629.78 


-1866.1  1510.0 

-1047.0  758.66 

1424.4  533.97 

53.925  2193.0 

63.402  -302.47 

-  93.805  656.54 
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We  want  to  determine  how  'close'  the  above  approximations  are  to  the  exact  solution. 
Thus  define  the  normalised  error  for  the  above  approximations  as 


£  = 


(3.10) 


where  A  is  the  exact  force  coefficient  matrix,/  is  a  load  vector,  and  E  is  the  error  matrix 
defined  as  E  =  A  -  A  where  A  is  some  approximation  to  A  (in  our  case  either  As%  or 
Aio%).  Finally  note  that  we  have  used  the  2-norm  in  the  above  definition,  where  the  2- 
norm  of  a  vector  x  is  defined  as  ||x||2  =  4x^x  .  The  above  normalised  error  measures  the 

accuracy  of  our  approximate  force  coefficient  matrix  A  compared  to  A  (the  exact  force 
coefficient  matrix).  We  can  now  investigate  the  effects  of  white  noise  on  error. 
Figure  3.1  shows  the  normalised  error  for  approximations  with  both  5%  and  10%  white 
noise  for  10,000  randomly  generated  unit  load  vectors  (these  random  vectors  were 
generated  as  discussed  earlier).  As  an  aside,  plots  of  the  normalised  error  distribution 
(akin  to  Figure  3.1)  using  100,  300,  and  1000  random  unit  load  vectors  (instead  of 
10,000)  show  excellent  agreement  with  the  distribution  resulting  from  the  10,000 
random  unit  load  vectors.  That  is,  the  distribution  appears  to  be  independent  of  the 
number  of  points  used  to  generate  it. 
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Figure  3.1  Sorted  normalised  error  versus  i th  random  vector  for  10,000  random 
vectors.  The  data  used  in  the  approximation  includes  either  5%  or  10%  white 
noise. 

As  we  would  expect,  Figure  3.1  shows  the  approximation  with  more  noise  produces 
larger  errors.  Figure  3.1  also  shows  that  the  regions  of  large  error  are  highly  localised 
since  only  a  small  portion  of  the  data  registers  large  errors.  For  example,  for  data  with 
10%  white  noise,  less  than  1%  of  the  data  have  error  exceeding  10%.  Consideration  of 
the  singular  values  shows  that  most  ill-conditioned  exact  matrices  A  would  yield 
localised  error  regions.  A  geometric  assessment  of  the  singular  values  also  leads  to  the 
conclusion  that  the  two  maximum  error  regions  become  more  localised  as  the 
dimensions  of  the  coefficient  matrix  increase,  especially  for  ill-conditioned  systems. 
(For  further  details  see  Appendices  E.l.)  Finally,  to  emphasise  the  point,  note  that  more 
than  95%  of  the  10% -white-noise  data  (98%of  the  5 %-whit e-noise  data)  has  less  than  5% 
error. 

Appendix  E  shows  that  the  normalised  error,  as  defined  in  Equation  (3.10),  for  any 
matrix  approximation  A  has  two  maxima  of  equal  magnitude  but  their  respective 
vectors  point  in  exactly  opposite  directions.  Using  this  result  we  can  safely  use  any 
maximisation  technique  to  produce  the  global  maximum.  Table  3.3  shows  this 
maximum  normalised  error  (and  associated  vector)  for  data  with  5%  white  noise.  The 
effect  of  the  number  of  data  points  in  developing  the  approximate  force  coefficient 
matrix  was  investigated.  The  first  column  shows  the  number  of  points  used  to  develop 
the  approximation  (or  in  other  words  the  approximation  order).  The  second  column 
shows  the  maximum  possible  error  for  any  unit  force  vector.  The  last  column  shows 
the  direction  of  the  unit  force  vector  corresponding  to  the  maximum  error.  As 
expected,  increasing  the  number  of  points  used  in  the  approximation  reduced  the  error, 
in  essence  averaging  out  the  white  noise.  Also,  contrary  to  expectations,  the  vectors 
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associated  with  the  maximum  error  do  not  appear  to  converge  with  an  increasing 
number  of  points  in  the  approximation.  This  small  sample  may  suggest  that  the  error 
matrix  E  is  uniformly  distributed  as  the  approximation  order  increases,  which  is  an 
assumption  we  made  in  Appendix  E.l  to  arrive  at  the  result  suggesting  that  the 
maximum  error  region  was  localised. 


Apx.  Pts. 

Max.  Error 

Direction  of  Max.  Error 

100 

3.386 

[-0.1452,  -9.357x10-4,  -0.1602,  0.1685,  0.2301,  -0.9337p 

300 

1.213 

[-0.6286,  -4.018X10-3,  -0.6897,  -0.04531,  0.08791,  -0.3454]T 

1000 

0.6686 

[  0.2604,  -3.522x10-3,  0.3524,  0.1478,  0.2134,  -0.8605]T 

3000 

0.2984 

[  0.1976,  8.507X10-4,  0.1858,  -0.1961,  -0.2249,  0.9151]T 

Table  3.3  Maximum  normalised  error  and  associated  vector  for  several  different 
approximations  using  a  different  number  of  points. 


A  graphical  representation  of  the  Table  3.1  results  is  shown  in  Figure  3.2.  This  figure 
shows  the  variation  of  normalised  error  (dotted  black  line)  versus  the  number  of  points 
used  in  the  approximation  of  the  force  coefficient  matrix  for  data  with  5%  noise. 


Number  of  Points  in  Approximation 

Figure  3.2  Plot  of  maximum  normalised  error  (dotted  black  line)  for  several 
approximations  of  the  force  coefficient  matrix  using  different  number  of  points  in 
the  approximation.  The  light  grey  line  shows  an  approximation  to  the  nonnalised 
error  upper  bound  (simply  given  by  the  ratio  of  singular  values). 

Again  we  see  that  the  maximum  normalised  error  decreases  as  the  number  of  points  in 
the  approximation  increases.  We  see,  however,  that  although  the  error  decreases  on 
average,  from  point  to  point  the  error  fluctuates  seemingly  randomly.  For  example, 
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using  a  240  point  approximation  the  maximum  normalised  error  was  found  to  be 
0.5238,  as  compared  to  an  error  of  0.6338  for  a  2617  point  approximation.  The  light  grey 
line  above  the  maximum  normalised  error  line  (dotted  black  line)  is  a  bound  obtained 
by  dividing  the  largest  singular  value  of  the  error  matrix  by  the  smallest  singular  value 
of  the  true  force  coefficient  matrix.  As  can  be  seen  this  maximum  bound  provides  a 
crude  approximation  to  the  maximum  error. 

Figure  3.3  shows  how  the  error  is  distributed  for  different  approximation  orders  using 
data  with  5%  white  noise.  The  10,000  random  points  were  sorted  by  the  magnitude  of 
the  normalised  error  so  as  to  obtain  the  distribution  plots.  Again,  the  error  decreases  as 
more  points  are  used  in  the  development  of  the  coefficient  matrix  approximation.  All 
approximations  show  good  agreement  with  the  exact  solution,  with  even  the  100  point 
approximation  registering  less  than  10%  error  for  more  than  98%  of  the  data. 


Figure  3.3  Normalised  error  for  10,000  random  unit  force  vectors  using  100,  300, 
1000,  and  3000  point  approximations  of  the  force  coefficient  matrix.  The  data  set 
used  to  develop  the  approximate  force  coefficient  matrix  contained  5%  white  noise. 


Although  we  have  excellent  approximations  of  the  coefficient  matrix,  remember  from 
Equation  (2.2)  that  we  need  to  invert  the  coefficient  matrix.  Performing  a  similar 
analysis  on  the  inverted  coefficient  matrix  produces  less  favourable  results.  Figure  3.4 
shows  the  normalised  error  distribution  (normalised  against  the  exact  inverse  of  the 
coefficient  matrix)  for  the  inverted  approximate  coefficient  matrices.  As  can  be  seen 
only  the  3000  point  approximation  gives  results  with  tolerable  errors.  This  result  might 
be  expected  from  consideration  of  the  unstable  procedure  of  matrix  inversion. 
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Figure  3.4  Normalised  error  for  1000  random  unit  force  vectors  using  the  inverse 
of  the  approximate  force  coefficient  matrices.  The  approximations  use  100,  300, 
1000,  and  3000  points  for  data  with  5%  white  noise,  and  1000  points  for  data  with 
10%  white  noise. 


Figure  3.5  Plot  of  maximum  normalised  error  (black  line)  for  several 
approximations  of  the  force  coefficient  matrix  using  different  number  of  points  in 
the  approximation.  Tire  dark  grey  line  shows  an  approximation  to  the  normalised 
error  upper  bound  (simply  given  by  the  ratio  of  singular  values).  The  light  grey 
line  shows  the  vector  norm  for  the  same  approximations. 
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Similar  to  Figure  3.2  we  can  generate  a  plot  of  maximum  normalised  error  versus  the 
number  of  points  in  the  approximation  for  the  inverse  approximate  coefficient  matrix; 
the  result  is  shown  in  Figure  3.5.  As  before  there  are  only  two  maxima  (both  the  same 
magnitude)  on  the  normalised  error  surface  generated  by  unit  load  vectors.  The  matrix 
maximum  was  thus  simply  found  by  maximising  this  error  surface.  The  normalised 
error  matrix  bound  is  again  the  ratio  of  the  largest  singular  value  of  the  error  matrix  to 
the  smallest  singular  value  of  the  exact  matrix. 

A  vector  norm  would  not  only  include  the  inverse  of  the  approximate  force  coefficient 
matrix,  but  would  also  incorporate  the  least  squares  approximate  of  the  stress  at  the 
zeroth  strain  gauge.  We  begin  by  defining  an  appropriate  measure  of  the 
approximation's  error.  From  Equation  (2.3)  the  stress  at  the  zeroth  gauge  was  given  by 

<t0  =  cqcq  +  a2<J2  h - \-aKoK,  which  may  be  written  in  vector  notation  as  cr0  =  aT  ■  a  , 

where  a  =  [cr,  a2  ...  aKJ  is  the  vector  of  coefficients  and  a  =  [cq  o2  ...  oKY  is 
the  vector  of  stresses,  for  the  gauges  1, 2 r...,K  respectively.  Given  an  approximation  d0 
to  the  exact  stress  a0  we  want  to  determine  some  sort  of  normalised  accuracy  of  that 
approximation.  We  cannot  simply  use  the  normalised  error  definition 

£max  =  max  (<J0  -  <70  )/<T0  since  the  denominator  becomes  zero  for  some  unit  stress 

IMM 

vectors  (given  by  ||cr|[  =  1 )-  However,  expanding  this  error  definition 
(aT -cr-a1 -cr)/(aT  cr)=[(d:-a)T  (j]/(aT  a),  where  a  is  the  cr0  approximation's 
vector  of  coefficients,  we  obtain  the  genesis  for  an  appropriate  normalised  error 
definition 


(a-a)T  (a -a) 


a  a 


(3.11) 


The  above  definition  may  be  thought  of  a  vector  norm.  Note  that  this  definition 
considers  only  the  coefficients,  and  is  thus  extremely  simple  to  calculate. 

Values  of  the  normalised  error  (vector  norm)  versus  the  approximation  order  are 
plotted  as  the  light  grey  line  in  Figure  3.5.  Note  the  close  relation  between  the 
maximum  error  of  the  normalised  matrix  (black  line)  and  the  vector  norm.  This 
suggests  that  most  of  the  error  comes  from  the  matrix  inversion  procedure  as 
suspected. 

Finally,  to  visually  gauge  the  accuracy  of  some  of  the  approximations  a  scatter  plot  (of 
exact  stress  versus  approximate  stress)  for  two  approximations  (1000  and  3000  point 
approximations)  is  shown  in  Figure  3.6.  This  figure  emphasises  the  result  shown  in 
Figure  3.4,  which  depicts  large  errors  for  the  1000  point  approximation. 


21 


DSTO-RR-0171 


5%  white  noise 


2000 ; 


Cfl 

W 

O) 

-M 

LD 

j 

1000  - 

O) 

"fS 

§ 

0 

X 

o 
<  ■ 

1 

a, 

< 

1 

-1000  - 

-2000 


I 


-2000  -1000  0  1000  2000 
Exact  Stress 


Figure  3.6  Scatter  plot  of  1000  random  unit  force  vectors;  the  horizontal  axis  shows 
the  exact  stress  value  while  the  vertical  axis  shows  the  approximated  stress  value. 
Approximations  of  the  force  coefficient  matrix  using  1000  (black)  and  3000  (grey) 
points  are  compared  for  data  with  5%  white  noise. 


3.3  Noisy  Stress  Estimations  with  Unknown  External  Loads  (vector 
technique) 

We  produce  some  unexpected  results  in  this  section,  namely  that  the  more  arduous 
method  of  the  previous  section  provides  better  results  than  stress  estimation  without 
the  force  coefficient  matrix. 

Figure  3.7  uses  the  normalised  error  definition  given  by  Equation  (3.11)  to  demonstrate 
the  relation  between  percentage  noise  in  the  data  and  normalised  error.  There  are  two 
disturbing  results  in  this  plot.  Firstly,  the  normalised  error  tapers  to  approximately 
unity  as  the  noise  increases  above  approximately  1%.  The  second  disturbing  result  is 
that  the  accuracy  of  the  approximation  appears  to  be  independent  of  the  number  of 
points  used  to  develop  that  approximation,  that  is,  accuracy  is  independent  of 
approximation  order! 
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Figure  3.7  Variation  of  normalised  error  with  percentage  noise  for  100,  1000,  and 
10,000  point  stress  approximations. 
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The  first  of  these  disturbing  results  may  actually  be  explained  on  closer  consideration 
of  what  is  happening  as  the  signal-to-noise  ratio  of  the  data  set  decreases.  As  we 
increase  the  noise  the  stresses  from  the  six  gauges  (used  to  approximate  the  stress  for 
the  zeroth  gauge)  become  uncorrelated  with  the  stress  at  the  zeroth  gauge.  Let's 
consider  a  simple  example  of  what  happens  to  the  linear  least  squares  fit  of 
uncorrelated  data. 


-0.5 


-1  -  : _ -L-i-i-'l-  ‘ 

-1  -0.5  0  0.5  1 

x 

Figure  3.8  Example  of  uncorrelated  data,  for  which  the  least  squares  linear  fit  is 
y=A+Bx.  As  the  number  of  approximation  points  tends  to  infinity  the  coefficients 
A  and  B  tend  to  zero. 
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Figure  3.8  shows  300  points  whose  x  and  y  coordinates  vary  randomly  (and 
independently)  between  -1  and  1.  The  best  linear  fit,  based  on  least  squares,  is  given  by 
y  =  -0.0277  -  0.0254x .  As  can  be  seen  (and  is  intuitively  obvious),  the  line  with  near 
zero  slope  best  fits  these  data.  If  the  amount  of  data  increases  and  the  x  and  y 
coordinates  becomes  more  uncorrelated,  then  in  the  limit  the  slope  approaches  zero. 
Referring  back  to  Equation  (3.11)  we  see  that  if  the  approximation's  coefficient  vector 
a  tends  to  zero  then  the  normalised  error  will  tend  to  unity.  This  fact  then  explains 
why  in  Figure  3.7  the  error  tends  to  unity. 

At  first  it  was  thought  that  the  original  linear  dependence  of  the  gauge  stresses  was  to 
blame  for  these  poor  results  as  soon  as  any  noise  entered  the  system.  Figure  3.9  plots 
the  normalised  errors  upon  deletion  of  different  gauges  compared  (as  a  ratio)  to  the 
normalised  errors  from  deleting  the  3rd  gauge.  Figure  3.9  verifies  the  fact  that  the 
smallest  normalised  errors  are  obtained  by  eliminating  the  3rd  gauge,  as  we  had 
predicted  using  the  condition  number  analysis  in  §  3.1. 


Figure  3.9  Comparison  of  errors  when  different  gauges  are  deleted.  The  comparison 
shows  the  ratio  of  the  error  from  deleting  the  nth  gauge  to  the  error  from  deleting 
the  3rd  gauge. 


Figure  3.10  compares  the  exact  solution  (the  horizontal  axis)  to  the  estimated  solution 
(the  vertical  axis)  using  1000  randomly  selected  stresses  (for  gauges  1-7  excluding  3) 
ranging  between  -10  and  10.  The  plots  on  the  left  and  right  show  the  results  when  the 
data  sets  used  to  develop  the  approximation  contained  5%  and  0.5%  white  noise 
respectively.  For  each  plot  two  approximations  were  developed,  one  using  1000  points 
(shown  as  black  dots)  and  the  other  using  3000  points  (shown  as  grey  dots).  Figure  3.10 
emphasises  the  fact  that  even  a  small  amount  of  noise  is  all  that's  needed  to  completely 
corrupt  the  approximation.  These  plots  also  show  that  the  stress  is  consistently 
underestimated  when  the  approximations  were  developed  with  data  containing  both 
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5%  and  0.5%  noise.  For  example,  for  the  0.5%  noise  plot  when  the  structure  is 
experiencing  2000  units  of  stress,  the  stress  estimated  is  less  than  1200  units  of  stress 
(on  average).  The  under-estimation  is  even  worse  for  the  5%  noise  plot.  Perhaps  this 
under-estimation  is  related  to  the  near  zero  coefficient  vector  of  the  approximation? 


Figure  3.10  Comparison  of  exact  solution  and  estimated  solution  for  1000  random 
data  points  with  5%  noise  (left)  and  0.5%  noise  (right),  using  1000  points  (black) 
and  3000  points  (grey)  approximations. 


In  order  to  gain  some  understanding  of  these  anomalous  results  we  investigate  a 
simpler  model  (yielding  a  2x2  matrix  system)  in  the  next  section. 
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4.  Preliminary  Results  of  a  Simplified  Truss 


In  this  section  we  develop  a  five-member  simple  truss  to  further  investigate  the 
properties  of  the  vector  and  matrix  solution  techniques. 

Figure  4.1  illustrates  the  statically  determinate  five-member  truss  we  will  use  in  this 
section.  All  joints  are  pinned,  joint  D  is  grounded  via  a  pin  restraint,  and  joint  B  is 
grounded  via  a  roller  bearing.  The  horizontal  distance  from  joint  D  to  joints  A,  B,  and  C 
are  1,  2,  and  1.1  respectively.  All  distances  have  been  non-dimensionalised  by  the 
vertical  distance  from  joint  D  to  joint  C,  and  thus  this  vertical  distance  is  unity.  We 
apply  a  horizontal  force  of  magnitude  Px  at  joint  C,  and  a  vertical  force  of  magnitude 
Py  at  joint  A.  We  will  determine  the  force  in  member  CD  using  the  forces  in  members 

AD  and  AB.  (The  choice  of  truss  and  forces  will  be  explained  in  a  future  report.  This 
section  merely  reports  preliminary  results  of  a  numerical  simulation.) 


Figure  4.1  Simplified  truss:  The  horizontal  distances  from  joint  D  to  joints  A,  B, 
and  C  are  1,  2,  and  1.1  respectively.  The  vertical  distance  from  joint  D  to  joint  C  is 
1. 


The  exact  forces  within  the  members  AD,  AB,  and  CD  are  given  respectively  by 


(4.1) 

(4.2) 
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and 


Fad  =  0.4500 Px  -  0.5500Py , 
Fab  =  0.4500P,,  -  0.4500Py , 
Fcd  =  0.7433P,  +  0.7433PV . 


(4.3) 
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Hence  we  can  solve  for  the  force  in  member  CD  in  terms  of  the  forces  in  members  AD 
and  AB  yielding  the  relation  FCD  =-14.8 7FAD  +16.52 FAB .  (Note  that  these  results  are 
quoted  to  four  significant  figures.)  The  above  relation  is  recovered  whenever  least 
squares  (LS)  is  used  to  fit  input  and  output  data  with  no  noise.  (The  input  data  are 
given  by  the  forces  FAB  and  FAD,  while  the  output  data  are  determined  from  the 
resulting  force  FCD.)  Using  noisy  data,  however,  we  find  that  as  we  increase  the 
number  of  points  in  the  LS  fit  the  vector  method  produces  large  errors.  On  the 
contrary,  the  matrix  method  results  improve  whenever  the  number  of  points  used  in 
the  LS  fit  is  increased. 

Figure  4.2  shows  the  distribution  (or  frequency  plots)  of  1000  LS  solutions  for  both  the 
vector  and  matrix  solution  techniques.  As  before  the  input  data  (loads  Px  and  Py )  were 

randomly  generated  using  a  uniform  distribution.  The  resulting  loads  in  members  AD, 
AB,  and  CD  were  determined  using  Equations  (4.1)-(4.3).  So  that  each  input  load  pair 
{Px,Py)  resulted  in  a  data  point  of  the  form  (FAD,FAB/FCD).  Each  LS  solution  used  128 

data  points  to  determine  the  fit.  Both  distribution  plots  divided  the  LS  solution  data 
into  16  data  bins,  and  on  a  logarithmic  scale,  all  bins  were  of  the  same  width.  In  a 
statistical  sense,  the  matrix  technique  produced  smaller  errors  than  the  vector 
technique  for  the  128  point  fit.  Both  the  vector  and  matrix  distributions  are  almost 
normal  on  a  logarithmic  scale.  The  slight  skewness  towards  zero  normalised  error  can 
be  explained  by  the  fact  that  the  error  is  bounded  from  above,  but  it  is  still  possible  to 
obtain  a  solution  with  zero  normalised  error. 
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Figure  4.2  Distribution  of  1000  least  squares  (LS)  solutions.  Each  LS  solution  used 
128  data  points  to  determine  the  fit. 

From  this  distribution  we  can  easily  determine  the  lowest  outlier,  first  quartile  ( Q1 ), 
the  median,  the  third  quartile  (Q3),  and  the  upper  most  outlier.  This  distribution 
information  allows  us  to  estimate  both  spread  and  skewness.  Repeating  the  above 
distribution  analysis  for  fits  using  2,  4,  8,. ..,  1024  point  approximations  we  can  produce 
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the  plot  shown  in  Figure  4.3,  which  shows  the  distribution  of  normalised  error  for 
several  approximation  orders.  Note  that  each  approximation  order  plotted  in  Figure  4.3 
consists  of  1000  normalised  errors  for  the  LS  fit. 
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Figure  4.3  Compariso7i  of  the  normalised  error  for  the  vector  and  matrix 
techniques.  Shown  are  the  distribution's  upper  and  lower  outliers  (outer  edges  of 
light  grey  region),  first  and  third  quartile  (outer  edges  of  dark  grey  region),  and 
median  (thick  line). 
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Although  the  spread  of  the  error  for  the  vector  technique  decreases  as  the  number  of 
data  points  increases,  the  error  converges  around  a  wrong  solution.  Thus,  in  a 
statistical  sense  (using  the  median),  the  error  actually  (slowly)  increases  as  the  number 
of  points  used  in  the  LS  fit  increases.  We  see  the  reverse  for  the  matrix  technique, 
where,  as  expected,  both  the  error  and  spread  decrease  with  increasing  information 
(again  in  a  statistical  sense). 

Preliminary  results  suggest  that  this  phenomenon  occurs  only  when  there  is  some 
collinearity  within  the  input  data,  and  is  (surprisingly)  independent  of  the  condition 
number  of  the  system.  Initial  results  also  suggest  this  drawback  in  using  the  vector 
technique  can  be  overcome.  The  modified  vector  technique  involves  using  the  lowest 
number  of  data  points  possible  in  the  LS  fit  (for  our  two-by-two  system  this  means  two 
data  points).  We  obtain  a  set  of  these  simple  LS  fits,  then  determine  the  median  of  this 
set.  Figure  4.4  plots  the  normalised  error  of  the  modified  vector  technique  for  several 
approximation  orders.  As  can  be  seen,  the  lowest  order  approximation  is  better  than 
the  high  order  approximations  by  almost  an  order  of  magnitude. 
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Figure  4.4  Plot  of  the  normalised  error  (for  the  modified  vector  technique)  versus 
least  squares  approximation  order.  Each  particular  approximation  order  point  in 
this  plot  represents  the  normalised  error  of  the  median  of  1000  least  squares 
solutions. 


Do  not  confuse  these  two  last  figures.  In  Figure  4.3  we  have  plotted  the  median  (along 
with  outliers  and  quartiles)  of  the  distribution  of  normalised  errors.  In  contrast.  Figure 
4.4  shows  the  normalised  error  of  the  median  of  the  distribution  of  LS  solutions  versus 
approximation  order.  This  modified  approach  also  works  for  the  matrix  method  and 
produces  better  results  than  the  original  matrix  procedure.  The  advantages  of  this 
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modified  procedure  are  that  it  is  faster  to  implement,  requires  less  information,  and 
produces  better  results.  A  proof  of  this  procedure  will  be  given  in  a  future  report. 
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5.  Discussion 


In  this  brief  section  we  summarise  and  expand  on  some  of  the  topics  already  discussed. 

The  non-linear  analysis  suggested  that  we  place  all  strain  gauges  away  from  buckling 
structures.  We  did  not  develop  general  guidelines  on  what  a  'safe'  distance  from  the 
buckling  member  would  be,  since  we  suspect  that  any  such  analysis  would  be  both 
complex  to  undertake  and  largely  problem-dependent.  As  such  the  choice  of  strain 
gauge  locations  requires  engineering  intuition  and  logic,  and  where  feasible  a  finite 
element  analysis  check  should  be  made. 

We  also  showed  that  non-uniqueness  effects  need  to  be  foremost  in  the  investigator's 
mind  when  deciding  which  locations  are  amenable  to  solution  by  the  stress  transfer 
function  method.  Essentially,  any  structure  for  which  we  have  more  unknown  loads 
than  load  equations  would  give  rise  to  a  non-uniqueness  problem. 

Linear  independence  of  strain  gauge  measurements  is  critical  if  a  stress  transfer 
function  is  to  be  developed.  However,  we  also  noted  that  linear  independence  was 
insufficient  when  noisy  measurements  were  involved.  We  need  to  further  investigate 
the  extent  of  system  ill-conditioning  that  can  be  tolerated  while  still  obtaining 
satisfactory  solutions.  That  is,  if  we  slightly  perturb  the  exact  force  coefficient  matrix 
(by  adding  noise),  is  it  likely  that  the  perturbed  matrix  will  be  close  to  an  ill- 
conditioned  matrix?  An  alternative  direction  of  research  might  be:  what  structural 
characteristics  cause  the  system  of  equations  to  become  ill-conditioned  in  the  first 
place?  We  have  already  mentioned  inappropriate  location  of  strain  gauges  as  one 
possible  suspect.  Are  there  others,  for  example  truss  geometry,  and  if  so  can  ill- 
conditioning  effects  be  quantified  or  at  least  ordered  (so  as  to  prioritise  effects)? 

An  initial  investigation  of  ill-conditioning  within  the  underlying  system  of  equations 
(that  is,  the  exact  matrix  A)  showed  that  the  matrix  technique  was  more  stable  than  the 
vector  technique.  For  well-conditioned  systems  both  methods  yielded  comparable 
results.  Preliminary  results  have  also  shown  that  a  modified  vector  technique  can 
produce  results  that  are  better  than  the  matrix  technique.  These  initial  results,  however, 
require  a  more  rigorous  generalisation  and  associated  validation,  but  go  some  way  to 
explaining  the  results  observed  for  the  vector  technique. 

One  reason  for  pursuing  the  vector  technique  over  the  matrix  technique  is  the 
simplified  coefficient  vector  construction  (no  knowledge  of  the  external  loading  is 
required).  Both  experimentally  and  in  the  final  implementation,  the  accurate 
application  and  measurement  of  external  loads  and  measurement  of  stresses  (the 
matrix  technique)  would  require  significant  additional  work  when  compared  to  taking 
only  stress  measurements  (the  vector  technique).  However,  the  highlighted  ill- 
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conditioning  problems  would  need  to  be  addressed  before  the  vector  technique  could 
be  successfully  implemented. 

We  know,  at  least  for  the  near  ill-conditioned  system  developed,  that  the  errors  from 
the  vector  technique  did  not  improve  when  the  approximation  order  was  increased.  It 
might  prove  profitable  to  pursue  ideas  based  on  the  bootstrap  method.  (For  a  concise 
overview  of  the  bootstrap  method  see  Press  et  al.  [35] .)  For  example,  would  combining 
several  low  order  approximations  improve  results,  or  could  the  data  be  reorganised  to 
produce  the  true  underlying  distribution? 

We  know  that  the  matrix  and  vector  techniques  are  based  on  the  same  fundamental 
solution  procedure  (least  squares),  and  that  the  final  forms  of  these  solutions  are  the 
same  (a  list  of  strain  gauge  coefficients).  Why  then  is  the  matrix  technique  more  stable 
than  the  vector  technique  (at  least  for  the  ill-conditioned  problem  in  this  report)?  The 
only  hypothesis  that  readily  comes  to  mind  is  that  the  greater  amount  of  information 
stored  in  the  development  of  the  force  coefficient  matrix  helps  to  stabilise  the 
subsequent  inversion  procedure.  However,  it  is  still  hard  to  imagine  that  under  some 
situations  the  indirect  matrix  technique  should  yield  results  that  are  superior  to  the 
direct  vector  technique. 

We  conclude  this  discussion  section  with  several  open-ended  questions  arising  from 
the  results  of  the  preceding  sections  and  the  discussion  within  this  section. 

•  How  does  noise  affect  the  system's  condition  number? 

•  How  do  structural  characteristics  affect  the  ill-conditioning? 

•  Why  is  the  error  independent  of  the  approximation  order  for  the  vector  technique? 

•  Why  did  the  vector  technique  under-estimate  the  results? 

•  Why  does  the  matrix  technique  appear  to  be  more  stable  than  the  vector  technique? 

•  Are  there  situations  under  which  the  vector  technique  gives  results  superior  to  the 
matrix  technique? 

•  Will  a  data  analysis  method  based  on  the  bootstrap  method  improve  results? 
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6.  Conclusion 

Consideration  of  non-linear  effects  showed  that  buckling  (whether  local,  as  in  beam 
flanges,  or  ordinary,  as  in  the  skin)  would  lead  to  a  bi-linear  stress  response  at  best.  As 
such  we  recommended  all  strain  gauges  set  up  to  determine  the  stress  response  of  any 
underlying  linear  (or  almost-linear)  system  should  be  located  away  from  the  load  paths 
of  buckling  structures.  We  were  also  able  to  demonstrate  that  indeterminate  systems 
still  remained  linear,  at  worst  complicating  linear  relationships  and  adding  constants  to 
the  solution  procedure.  We  also  showed  that  the  under  some  conditions  a  static  stress 
model  would  yield  non-unique  solutions  if  measurements  were  made  outside  the 
relevant  structure.  For  a  simply  supported  beam,  we  also  determined  when  the  static 
stress  distribution  model  would  provide  a  poor  approximation  of  the  dynamic  stress 
distribution. 

We  examined  two  solution  techniques,  the  matrix  and  vector  techniques.  We 
developed  the  stress  transfer  function  (STF)  for  the  matrix  and  vector  techniques 
respectively,  with  and  without  external  load  information.  Although  the  matrix 
technique  was  susceptible  to  noise,  the  error  in  the  STF  approximation  could  be 
reduced  by  increasing  the  approximation  order  (that  is,  the  number  of  load-stress 
points  in  the  least  squares  fit).  In  contrast,  increasing  the  approximation  order  for  the 
vector  technique  had  an  insignificant  effect.  For  noisy  data,  the  vector  technique's  STF 
improved  only  when  the  signal-to-noise  ratio  was  well  above  realistic  levels.  Initial 
results  seem  to  suggest  that  the  poor  performance  of  the  vector  technique  was  related 
to  the  ill-conditioned  system  of  equations  governing  stress  in  the  truss.  The  matrix 
method,  however,  seemed  far  less  susceptible  to  this  ill-conditioning.  Finally, 
preliminary  results  also  suggest  that  taking  the  median  of  a  large  number  of  low  order 
approximations  yields  good  results  for  both  the  matrix  and  vector  techniques. 
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A.  Statically  Indeterminate  Structures 
Remain  Linear 


In  this  appendix  we  demonstrate  two  sets  of  examples  showing  that  statically 
indeterminate  structures  have  additional  complexity  (compatibility  equations  have  to 
be  additionally  solved),  but  the  problem  remains  linear.  In  fact  the  system  will  remain 
linear  unless  it  experiences  buckling,  large  strains  (that  is,  geometric  non-linearity),  or 
the  material  behaviour  is  non-linear  [36,  37], 

The  two  conditions  shown  in  Figure  A.l  demonstrate  how  redundantly  constrained 
structures  may  be  thought  of  as  determinate  structures  that  are  pre-stressed.  Thus  for 
stress  relationships  between  different  members  of  a  truss  system,  these  types  of 
redundancies  merely  change  the  constants  in  the  linear  system  of  equations. 


Figure  A.l.  Some  redundant  structures  (shown  on  the  left)  may  be  thought  of  as 
determinate  but  pre-stressed  structures  (shown  on  the  right). 


Figure  A.2  shows  a  truss  system  with  a  redundant  member  and  a  vertical  and 
horizontal  force  at  the  joint  labelled  A  (all  joints  are  pinned).  Using  indeterminate  truss 
analysis  the  forces  in  the  members  AD,  AB,  and  AC  are  given  respectively  by 


and 


1  I 

F*  +  V 2Rac/ 

R-ab  ~  +  2  Fy  ~  2^/2  Fac  ' 
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where  l/cAC  is  equal  to  the  product  of  area  and  Young's  modulus  for  the  member  AC, 
that  is  cAC  =1/(AacEac)  (a  similar  notation  holds  for  the  other  members).  The  last 
equation  was  derived  using  the  compatibility  equation  eAC  =  (eAB  +  5ead)/A,  where 
eac  is  the  strain  of  the  member  AC.  This  simple  example  demonstrates  that  the  stresses 
in  these  members  are  still  linear  functions  of  the  external  forces  Fx  and  Fy,  and  thus 
there  is  a  linear  relation  between  stresses  in  the  truss  members.  So  that  if  we  assume 
the  same  Young's  modulus  for  all  three  member,  then  the  tensile  stress  in  member  AC, 
for  example,  in  terms  of  the  stresses  in  members  AB  and  AD  is 

°  AC  =~^iaAB  +5(7  AD), 

where  Gab,  Gac,  Gad,  Aab,  Aac,  and  Aad  are  the  tensile  stresses  and  the  areas  in  members 
AB,  AC,  and  AD  respectively. 


1 


Figure  A.2.  Even  with  a  redundant  member  the  truss  system  is  still  linear. 
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B.  Solution  of  a  2-D  Statically  Determinate 
Helicopter  Truss 


In  this  appendix  we  solve  the  linear  system  of  equations  governing  the  forces  and 
moments  in  the  simple  two-dimensional  model  of  a  helicopter  truss  shown  in 
Figure  2.1  (re-drawn  below  in  Figure  B.l  with  additional  joint  labelling  and  the 
accelerations  shown). 


Unless  otherwise  stated  the  positive  sense  of  forces  and  moments  are  in  the  directions 
specified  by  the  main  and  tail  rotor  loads.  Members  AC,  BC,  and  CH  are  pinned  to  the 
slide  bearing  which  passes  through  member  DE.  Members  AB  and  DE  are  single  beams 
(that  is,  there  are  no  hinges  along  these  beams).  As  stated  earlier,  joint  A  is  pinned  to 
the  ground  and  joint  B  is  connected  to  a  roller,  while  tuac  and  hiad  are  point  masses. 
The  location  of  a  joint  or  strain  gauge  is  given  in  terms  of  its  x  and  y  coordinates 
subscripted  with  the  joint  letter  or  strain  gauge  number.  For  example,  X2  is  the  x- 
coordinate  of  strain  gauge  G2.  Angles  are  denoted  by  0  with  a  three  letter  subscript 
which  corresponds  to  joints  that  enclose  the  angle.  For  example,  0dac  is  the  angle 
formed  by  the  members  DA  and  AC  at  joint  A. 

The  solution  of  the  static  truss  is  shown  below,  where  all  reaction  forces  are  given  in 
the  positive  direction  (shown  by  the  axes  at  the  origin,  joint  A),  moments  are  positive  in 
the  anti-clockwise  direction,  and  all  member  forces  are  given  as  tensile  forces. 
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B.l  Main  Rotor  Loads  on  Cabin  Section 

In  this  subsection  we  determine  forces  and  moments  (and  some  shearing  forces)  on  the 
cabin  section  of  the  helicopter  (beams  AB,  AC,  BC,  and  DE )  from  main  rotor  loads  only 
(that  is,  Px/  Py,  and  Pm). 

The  constraint  reactions  from  the  ground  on  pin  A  and  on  roller  B  are  respectively 

Rax=~Px,  (B.1) 

(B.2) 

Rliy  ~  {Px  J/e  ~PyXD  ~  Pm)/XB  •  (B-3) 

Separating  member  DE  from  the  truss  and  considering  its  constraining  forces  gives 

RCx  =(Pm~  Px  Ve  )/yc  > 

RDx  =  ~PX  ~  RCx  ' 


B*Ay  Bgy  Py  , 


and 


and 


RDy  -  By  ■ 


In  addition  to  tensile  or  compressive  forces,  beams  AB  and  DE  (and  hence  gauges  Go, 
Gz,  and  G4 )  also  have  bending  moments  and  shear  forces.  (Shear  forces  at  the  strain 
gauges  are  not  shown  below  since  these  forces  do  not  influence  gauge  readings.)  The 
two  angles  0D ac  and  0dbc  may  be  written  down  in  terms  of  lengths  yielding 

6DAC  =  tan-1  (yc jxD )  and  6mc  =  tan"1  \yc / (xB  -xD)]. 

Then  the  forces  in  members  AB,  AC,  BC,  and  DE  are  respectively 

PAC  =-Fbc  sin0DBC/sin0DAC  , 

RBC  =  P-Cx  S^n  ®DAC  /Sin  (@DAC  +  ^DBc)' 

RADx  ~  ~RAC  COS  dD AC  ~  P-Ax  ' 

F ADy  =  -FaC  S^n  ®DAC  ~  RAy  ’ 


Fbdx  ~  Fbc  cos  6 , 


DSC  ' 


and 
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FBDy  ~  R-B)/  +  ^BC  S^n  ®DBC  ' 

where  the  forces  Fadv  and  FBnv  are  shear  forces  (positive  up). 

The  forces  at  the  gauges  may  be  similarly  calculated  to  give  tensile  forces  (subscript  T) 
and  bending  moments  (subscript  M),  where  the  subscript  number  denotes  the 
corresponding  strain  gauge.  So,  for  example,  Ft4  represents  the  tensile  force  at  gauge 
G4,  while  Fmo  denotes  the  (anti-clockwise)  bending  moment  at  gauge  Go.  The  forces  and 
moments  at  the  strain  gauges,  due  to  the  main  rotor  loads,  are  then 


O 

II 

Fmo~Fx(}/e  y0)+RaCVc  3/o)  Fm, 

(BA) 

Fn=FAo 

Fmi  =  ^  > 

(B.5) 

Ft2  =  FgDx  ' 

Fm2  =  ~FBD}/  (xb  ~X2J' 

(B.6) 

Ft3  ~  FBC> 

Fm3  =  0  / 

(B.7) 

Ft  4  =  FADx  , 

Fm4  ~  ~FaD;/X4  ' 

(B.8) 

0 

ll 

in 

u!r 

Fm  5  =  0. 

(B.9) 

B.2  Accelerations  on  Cabin  Section 

The  loads  and  moments  on  the  cabin  section  arising  from  mass  accelerations  are  given 
in  this  section.  The  ground  constraints  resulting  from  the  acceleration  of  mass  mAo  are 
denoted  by  hatted  variables  (A),  and  similarly  constraints  resulting  from  mass  mAc  are 
denoted  by  tilded  variables  (~).  The  mAo  and  mAc  resulting  reaction  forces  are  given 
respectively  by 


and 


B-Ax  —  mAD8x ' 

RAy  =~F-By  ~  mAD&y  ' 

Bgy  =  ~mAD&y  ( XMad  /Xb)> 

B Ax  ~  ~~^AC  8x  ' 

R-Ay  ~  -RBi/  ~  mAc8y  ' 
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RBy  =  mAc  igx  sindDAC  -  gy  COS0DAC 
where  xM«d  is  the  x-coordinate  of  the  point  mass  mAo.  Member  forces  easily  follow 

giving 

Fab  =  ~Fbc  cos6dbc  , 

and 

FBC  =  -  RBy  /sin  &DBC  ■ 

The  strain  gauge  forces  due  the  combined  masses,  which  we  denoted  by  hatted  tilde 
(£)  variables,  are 


F  =0 

rT  0  u  ' 

Fmo  _  ^ ' 

(B.10) 

Fji  =  Fbc  cos(Gdac  +  9dbc  ) , 

i~Mi  -  Fbc  sin(0DAC  +  0DBC )  , 

^COS0DAC  J 

(B.ll) 

Ft2  ~  Fab  , 

Fm2  ~  ~FBy  (xB  -  x2  ), 

(B.12) 

Ft  3  =  Fbc  > 

Fm3  =  0  ' 

(B.13) 

Ft 4  =  FAb  ' 

Fm  4  "  ~^By  iXB  ~  X4  )  / 

(B.14) 

o 

ll 

in 

<?u!r 

i*M5  =  0  ’ 

(B.15) 

Tail  Loads 

All  reactions  to  tail  loads  have  a  dashed  (')  notation.  Considering  tail  rotor  loads  only, 
the  constraint  reactions  from  the  ground  on  pin  A  and  on  roller  B  are  respectively 

R'  -  -T 

1KAx  —  1x' 

R'  -  -T  -  R' 

IMi/  ~  1  y  lxBy  ' 

and 

FBy  =  -Ty  +  Fch  [vc  cos0CHXX  (xB  xc)sin0  CHXX  ]/XB  ' 
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where  0CHXX  =  tan  1  [(yc  -yH)/{xH  - xc )]  is  the  angle  made  by  the  member  CH  with  the 
horizontal  (denoted  by  XX)  and  F'ch  is  defined  below.  To  help  with  the  derivation  of 
member  forces  Fbg  and  Fbh  below  note  that  the  geometric  ratios  yd  xb  and  ( xb-xc )/ xb  in 
the  expression  for  R'By  shown  above  have  the  equivalent  trigonometric  forms 

yc/xB=  tan dDAC  tan 0DBC /(tan dDAC  +  tan 0DBC ) 

and 

{xB-xc)/xB  =tan  6DAC /(tan  6 

DAC  +  tan0DBC). 

The  member  forces  are 

TxyF~Tv(xF-xB)-Tm 
yH  cos  9chxx  +(xh  —  xb)  sin  0CHxx 

Fch(cos@chxx  tan dDBC  -  sinOCHXX) 
sin  0DAC  +  cos  dDAC  tan  6mc 

,  _  -Fch  cosdCHXX  +Fac  cos0DAC 

BC  —  _ r\  ' 

cos  6DBc 

„  _  Ty  -  Tx  tan  0BHXX  +  FCH  cosdCHXX  (tan  0BHXX  +  tanflCHXX ) 
cos  0/jgxx  (tan  @bgxx  ~  tan  9Bijxx  ) 

_Ty—Tx  tan  0  BCXX  +  FCH  cos0cwvx(tan  Q  BGX X  +  tan  0CHXX ) 
cos  0BHXX  (tan  6  BHXX  —  tan  8  BGXX  ) 


f~AB  ~  F AC  COS  ®DAC  ^Ax  ■ 

And  thus  the  gauge  moments  and  tensile  forces  are 


f'  -  n 

rT  0  —  u  ' 

F'  -  0 

rM  0 

(B.16) 

F'  -  F' 

rT 1  _  rAC  ' 

& 

>-» 

II 

O 

(B.17) 

F'  -  F' 

r T2  ~  rAB ' 

(B.18) 

F'  -  F' 

rT3  rBC ' 

^M3  ~  0  / 

(B.19) 

F'  -  F' 

rT 4  rAB  ' 

^M4  =  0  ' 

(B.20) 
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and 


f'  -  V' 
rT  5  -  rCH  ' 

O 

II 

in 

h? 

(B.21) 

f'  -  C' 
f  T6  ~  rBH  ' 

^M6  ~  0  ' 

(B.22) 

F'  -  F' 

rT?  — r BG, 

Fm7=0. 

(B.23) 

B.4  Stresses  at  Strain  Gauge  Locations 


Now  that  we  have  all  the  member  tensile  forces  and  moments  under  the  loading 
shown  in  Figure  B.l  we  can  simple  state  the  stresses  at  the  strain  gauge  locations. 
Tensile  forces  F  will  cause  a  tensile  stress  a  given  by  the  relation  cr  =  F/A ,  where  A  is 
the  beam's  cross  sectional  area  at  the  strain  gauge  location.  Similarly  the  bending 
moment  Fm  will  cause  a  stress  o  given  by  a  -  FMr/I ,  where  I  is  the  beam's  second 
moment  of  area  at  the  strain  gauge  location  and  r  is  the  distance  from  the  neutral  axis 
to  the  gauge-beam  contact  surface.  We  then  have  in  general  that  the  stress  at  the  ;th 


strain  gauge  is 


a’  A 


1  ( 


7  V 


+  F;  +  F;  + 


F Mj  +  +  ^Mj 


\ 


where  Aj  is  the  cross  sectional  area,  r,  is  the  distance  from  the  neutral  axis  to  the  strain 
gauge,  and  Ij  is  the  second  moment  of  area  at  the  ;th  strain  gauge  location.  (The  sign  in 

front  of  the  moments  FMj  ,  FM; ,  and  F'Mj  is  positive  if  the  ;th  strain  gauge  is  in  tension 
under  that  moment  and  negative  otherwise.)  More  concretely  we  have  that 

o0=~-+^[px(yE-yo)+Rcx(yc-yo)-pm]' 

*0 


o1  =— [f 


AC  +  ffiC 


COS  {dDAC  +  ®DBC 


)+F'ac ] 


+  4LF 


h 


BC 


Xr  ~X, 


cos  6 


D  AC 


sin(0DAC  +  0DBC 


l 


c2  - UCBDx  +  ^AB  +  )+  r2  (  ^BDy  X2  )  ' 

A2  i2 

03  =  “7  (^BC  +  ^BC  +  F BC  )' 

A3 

04  ~  ~2  (^ADx  +  ^AB  +  ^AB  )+  ^ADyX4  —  ^By  (XB  —  X4  )]' 

A4  i4 
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^5 


and 


a7 


where  beam  forces  and  ground  reactions  are  given  by  Equations  (B.1)-(B.23). 

The  stress  equations  shown  above  were  verified  numerically  using  NASTRAN,  the 
results  are  shown  in  Table  B.l.  The  last  column  gives  a  numerical  comparison  between 
the  analytic  results  and  the  NASTRAN  model  in  terms  of  relative  error,  defined  as 
(NASTRAN-analytic)/ analytic.  The  input  loading  for  this  comparison  was  arbitrarily 
chosen  to  be  Px=7,  Pv= 8,  and  Pm- 9  for  the  main  rotor,  Tx= 3,  T„=l,  and  Tm=13  for  the  tail 
rotor,  and  gx= 5  and  gy= 6  for  the  accelerations  (or  inertial  loads  as  termed  in 
NASTRAN).  As  can  be  seen  the  derived  stresses  compare  well  with  the  NASTRAN 
model.  For  some  reason  the  stress  recovery  at  strain  gauge  Go  gave  the  incorrect  result 
-6.7528xl06.  However,  the  correct  result  is  shown  in  Table  B.l  derived  from 
NASTRAN's  beam  axial  force  and  bending  moment  at  strain  gauge  Go.  Increasing  the 
number  of  elements  by  a  factor  of  four  (by  twice  splitting  each  element  in  half) 
eliminated  this  NASTRAN  self-inconsistency  between  the  recovered  stress  and  the 
axial  force  and  bending  moment  at  gauge  Go. 


Strain  Gauge 

Analytic 

NASTRAN 

Relative  Error 

0 

-1.2725x106 

-1.2725xl06t 

-2.3477x10-6 

1 

3.6824x106 

3.6824x106 

3.5439x10-8 

2 

2.9202x106 

2.9202x106 

1.7874x10-6 

3 

3.9764xl04 

3.9771x104 

-1.7304x10-4 

4 

6.8519x106 

6.8519x106 

1.2394x10-6 

5 

-4.1779x104 

-4.1780x104 

-1.8601x10-5 

6 

-7.8983x104 

-7.8983x104 

-3.6965x10-7 

7 

1.0504x105 

1.0504x105 

-4.1463x10-7 

Table  B.l  Comparison  between  the  analytic  results  and  a  NASTRAN  model. 


Joint 

A 

B 

C 

D 

E 

F 

G 

H 

Point 

i 

2 

3 

4 

5 

6 

7 

8 

X 

0.0 

3.0 

2.0 

2.0 

2.0 

7.5 

3.6 

3.3 

y 

0.0 

0.0 

1.0 

0.0 

1.5 

0.8 

0.6 

0.9 

Table  B.l  Location  of  NASTRAN  truss  model  joints. 


t  This  value  was  derived  using  the  NASTRAN's  beam  axial  force  and  bending  moment  instead 
of  the  NASTRAN's  beam  stress  recovery  option,  which  was  used  for  all  the  other  strain  gauges. 
See  main  text  for  further  details. 
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Line 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

Point  1 

1 

1 

4 

2 

2 

2 

8 

8 

8 

7 

Point  2 

2 

3 

5 

3 

7 

8 

3 

7 

6 

6 

Num.  of  Elements 

6 

4 

3 

3 

2 

2 

3 

1 

8 

7 

Table  B.3  Number  of  elements  in  each  beam  for  the  NASTRAN  truss  model. 


The  NASTRAN  model  used  beam  (BAR2)  elements  that  were  pinned  using  multi-point 
constraints.  The  x-  and  y-coordinates  of  the  truss'  joints  are  shown  in  Table  B.2,  while 
Table  B.3  shows  the  number  of  elements  in  each  beam.  All  beams  had  the  same, 
circular  tube  cross-section  (of  inner  radius  0.015  and  thickness  0.005),  Poisson  ratio 
(v=0.3),  and  Young's  modulus  (E=200xl09).  The  two  point  masses  were  wac=9  and 
mAD=7. 
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C.  Effect  of  Vibration  on  Static  Assumption 


In  this  appendix  we  determine  under  what  conditions  the  assumption  of  negligible 
vibration  leads  to  grossly  inaccurate  results.  We  begin  by  investigating  the  dynamic 
properties  of  a  simple  beam  pinned  at  both  ends,  as  shown  in  Figure  C.l,  under  a 
transverse  dynamic  loading. 


Figure  C.l.  Transverse  vibration  of  a  pinned  beam  with  constant  mass  and 

stiffness. 


The  beam  shown  above  is  governed  by  the  following  partial  differential  equation 
[38,  p.  222] 

mpt  +  EIpt  =  f(x,t),  0<  x<L, 

dt~  dx 


where  y(x,  t )  is  the  vertical  displacement  at  time  t  and  at  a  point  x  along  the  beam,  m 
and  El  are  the  mass  and  stiffness  (both  constant)  of  the  beam,  and  f(x,t )  is  the  loading 
distribution  as  a  function  of  space  and  time.  Consider  the  initial  conditions 


y(x, 0)  =  y0  (x)  (initial  displacement). 


dt  t=0 


v0(x) 


(initial  velocity). 


and  boundary  conditions 


y(  0,0=0, 


y(L,t)=0/ 


El 


El 


?y 

dx2 

z2y 


=  o. 


x=0 


dx2 


=  0, 


x=L 
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where  the  second  condition  on  both  lines  is  equivalent  to  a  zero  bending  moment  since 
M  =  EId2y/dx2 .  The  solution  for  transverse  vibration  of  the  beam  shown  in  Figure  C.l 
under  the  above  initial  and  boundary  conditions  is  given  by  Meirovitch  [38,  pp.  235- 
37]  as 

y(x,  t)  =  Y  Y„  (x)|  —  [  Q„  (r)sin[fl)„  (t  -tjfrr  +  q  cos(cont )+ sin(ft)„f)l, 

ti  [®n  J0  J 

where  the  natural  frequencies  and  modes  are  given  respectively  by 


and 


and  the  generalised  forces,  coordinates,  and  velocities  are  given  respectively  by 


and 


Qn(t)  =  \  f(X’t)Y„(x)dx , 

o 

L 

q„B  =  jmy0(x)Yn(x)dx , 
o 

L 

i ln0  =  \mv0(x)Yn(x)dx . 
0 


We  want  to  determine  how  the  stress  response  of  a  beam  varies  as  the  forcing 
frequency  approaches  the  first  natural  frequency  of  the  beam.  Under  a  static  force  F 
applied  at  a  point  x  =  aL  along  the  beam,  where  as  (0,1),  the  static  stress  is  given  by 
Young  [37,  p.  100]  as 

<T9tat  =  -r-  [(1  -  a)x  -(x-  aL)h(x  -  aL)], 

where  z  is  the  distance  from  the  neutral  axis  to  the  fibre  at  which  stress  is  to  be 
measured  (usually  the  extreme  fibre),  x  is  measured  from  the  left  end  of  the  beam,  I  is 
the  moment  of  inertia  with  respect  to  the  neutral  axis,  and  h(x  —  aL)  is  a  step  function 
such  that  h(x  -  aL)  =  0  if  x  <  aL  and  h(x  -  aL)  =  1  if  x>aL. 
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We  want  to  compare  the  above  static  solution  with  that  of  a  cyclically  applied  load  and 
corresponding  vibratory  response.  Assume  the  distributed  loading  shown  in 
Figure  C.l  is  given  by 


f(x,  t )  =  -F  8(x  -  ah)  sin (coFt  +  dF ), 


where  8(x-aL)  is  the  delta  function  defined  such  that  f  8(x-r)f(x)dx  =  f(z)  for 

J  z-E 

any  e  >  0 .  Thus  the  above  forcing  function  represents  a  cyclic  point  force  at  x  =  aL 
having  a  maximum  force  F,  angular  frequency  £Uf ,  and  phase  shift  6r  . 

For  simplicity  we  consider  the  case  of  zero  initial  conditions,  y0(x )  =  v0(x)  =  0,  and  a 
zero  phase  shift  dF  =  0  which  yields  the  solution 


where 


y(x,t)=- 


F  (  2L3 
I  7T4E 


sin (n ;tqi  )  sin (n n  xj L ) ^ 
n4  -  A2 


=  sin(tuFf)--^-sin| 


fn2 


and  /i,  =  coF/a>:  is  the  ratio  of  the  forcing  frequency  to  the  first  natural  frequency.  So 
that  as  /?,  tends  to  zero  the  forcing  is  slow  when  compared  to  the  unforced  vibration  of 
the  beam,  while  when  /3,  tends  to  unity  the  forcing  frequency  is  close  to  the  first 
natural  frequency. 

The  relation  between  stress  and  bending  moment  is  o  -  Mz/I ,  while  M  =  EId2y/dx2 
gives  the  relation  between  bending  moment  and  displacement.  Thus  twice 
differentiating  the  displacement,  given  above,  yields  the  dynamic  stress 

_Fzf2L ^  (pn(t )  sin(nna)sin(nK  x/L) 

^ =  q?  B?[i  -fa/sj]  ?  ' 


The  series  above  can  easily  be  shown  to  be  convergent.  Clearly  the  expression  after  the 
summation  symbol  in  the  above  series  is  less  than  2/n2  ,  in  modulus,  for  n  sufficiently 

large.  And  since  the  ^jT  2 jn2  is  absolutely  convergent,  then  the  series  in  the  expression 

of  dynamic  stress  is  also  absolutely  convergent.  Thus  we  may  re-arrange  the 
summation  of  the  series  as  follows. 


We  expect  the  series  to  be  dominated  by  the  first  term,  so  we  separate  the  first  term  out 
giving 
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sin  bra)  sin  (ttx/L)  1  sin  (n^a)  sin 


-A 

Using  Prudnikov  [39,  p.  743  #4]  we  know  the  solution  of  a  series  similar  to  the  first 
series 


sin  sin  {coFtn2/PiJ  sm(nxa)sin{njc  x/ L) 

“FftT  (  '  (  v  }  6  t-fc/.1)1] 


(nnxjL) 


*^dyn 


l  sin(&)f  f) 

K*  ) 

Y  s^n(n7ra)sl^n(n7r -  -sin(^a)sin(^ x/l)+ — [(1  - ct)x  - (x - al)h(x  - aL)] , 

U  n  21 

and  we  can  approximately  bound  the  second  series  by 
^  sjn (coFtn2/ fij)  sm(nna)sin(nK  x/L ) 

S  fl  »* 


< 


y  1  _f4 

^“90 


1. 


This  is  only  an  approximate  bound  because  the  term  1  -  (/J3/n2)~  is  less  than  unity.  But 
since  P:  <  1  this  term  rapidly  approaches  unity  as  n  increases,  and  hence  the  above 
bound  is  an  excellent  approximation  to  the  true  bound.  The  same  argument  applies  to 
the  first  series.  In  fact,  the  modulus  of  the  series  shown  above  has  the  exact  bound 
(7T4/90-l)/|l-A/16|. 

Using  the  above  two  relations  we  may  approximate  the  dynamic  stress  by 


<xdyn  =y^j|siriK0  sin(7ra)sin(;r  x/l|^— - 1  +~  [(1  ~a)x  ~  (x  ~  aL )  h(x  ~  aL^ 


sin 


(i-ft2) 


sin(;ra)sin(7r  x/L)+  c1 


(  re4 
90 


-1 


where  c3  e  [-1,1] .  Rewriting  the  above  expression  in  terms  of  the  static  stress  and 
setting  c3  to  zero  (see  discussion  below)  we  obtain 


where 


and 


Cdyn  =  CTstat  |sm(fflFf)  +  gl  («/  X/L)  g2  (Pi  > 

a  (n  r,r\_(  2  "I  sin(^)sin(^x/L) 

^  '  l^2  J[(l-a)x/L-(x/L-a)h(x/L-a)] 
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g2(P1,(oFt)=  Pl-2  \fi,  sin(wff)-sin(ft)F  !/&)]. 

l  Pi 

Surface  plots  of  the  functions  g1  and  g2  are  shown  in  Figure  C.2  and  Figure  C.3 
respectively.  The  plot  of  function  g{  shows  that  the  dynamic  stress  and  static  stress  are 
similar  when  the  excitation  and  reading  locations  are  close  together  (that  is,  x~aL) 
and  most  similar  when  the  excitation  and  the  reading  location  are  at  either  end  of  the 
beam. 


Figure  C.2.  Function  g1  shows  the  dynamic  stress  is  very  close  to  the  static  stress 
when  the  excitation  and  reading  location  are  close  together. 

Finally  by  taking  the  limit  of  g^e)  and  gx(e,  1  -  e)  as  e  tends  to  zero  yields  the  lower 
and  upper  bounds  on  function  gx  respectively 

0<gi<2. 

As  expected  the  plot  of  the  function  g2  shows  that  the  dynamic  stress  is  closest  to  the 
static  stress  when  /},  is  small  (that  is  when  the  loading  frequency  is  much  smaller  than 
the  first  natural  frequency).  It  is  easy  to  show  that  the  function  g2  is  bound  by 
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Pi 

i-A 


—  82  — 


_A_ 

i-A' 


thus  we  see  that  resonance  results  as  /3,  approaches  unity  (again  as  expected).  The 
other  effect  to  notice  in  the  plot  of  g2  is  the  chirp  effect  as  coFt  grows  (this  chirping 
effect  is  shown  more  clearly  in  Figure  C.4). 


Figure  C.3  Function  g2  shows  chirping  in  the  /3,  direction  as  coFt  grows. 


Figure  C.4.  Chirp  effect  of  the  function  g2  evaluated  at  coFt  =  4k  . 
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Numerical  experiments  verify  that  assuming  c1  =  0  is  a  good  approximation.  To  see 
why  let's  examine  the  original  series  that  contained  the  term  sin (coFtn2/ fa).  For  n 
large  we  see  that  this  trigonometric  function  varies  rapidly  between  [-1,1]/  and  hence 
acts  almost  as  a  random  number  between  these  ranges.  Thus  we  would  expect  that  for 
a  large  number  of  cases  this  would  bring  the  original  series  close  to  zero.  Setting  c1  to  1 
and  -1  results  in  an  approximate  bound  (which  may  be  thought  of  as  an  error  bound) 
of  the  exact  solution's  approximation. 
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Figure  C.5.  Comparison  of  the  'exact'  and  approximate  scaled  dynamic  stresses  for 
500  points  having  randomly  distributed  arguments.  The  bounds  are  determined  by 

setting  c1  to  -1  and  1. 


Figure  C.5  shows  the  'exact'  solution,  calculated  by  summing  over  the  first  100  terms  of 
the  series,  and  the  upper  and  lower  bounds  of  the  approximate  solution.  The  500  points 
were  obtained  by  setting  coF  =  k  and  L  =  1  and  independently  randomly  distributing 
t ,  a ,  and  x  between  [0,1].  For  clarity  the  resulting  set  of  points  were  sorted  by  the 
magnitude  of  the  'exact'  solution  and  plotted  (this  explains  the  smooth  monotonic 
growth  of  the  'exact'  solution  of  a  randomly  selected  set).  As  can  be  seen  the  upper  and 
lower  approximate  bounds  are  extremely  accurate.  Also  the  'exact'  solution  appears  to 
fall  at  the  mid-point  of  the  upper  and  lower  bounds,  justifying  the  approximation 
c1  -  0 .  There  was  some  concern  about  the  accuracy  of  the  approximation  as  /?, 
approached  unity,  since  the  series  approximations  were  made  assuming  /3,<§C  1 .  A  plot 
of  the  scaled  dynamic  stress  with  all  variables  set  to  the  same  values  as  for  Figure  C.5, 
except  for  /?,  =  0.95 ,  showed  a  similar  accuracy  in  the  upper  and  lower  bounds  for  the 
approximation  of  dynamic  stress  as  shown  in  Figure  C.5. 
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Using  the  two  bounds  on  the  functions  gy  and  g2  we  arrive  at  an  approximate  bound 
on  the  difference  between  the  dynamic  stress  and  the  static  stress  as  a  function  of  the 
forcing  frequency  to  the  first  natural  frequency  ratio 


<Uvn 

^stat 


< 


1  + 


2A 

i-A  ‘ 


However,  if  we  choose  a  appropriately,  then  the  bound  on  function  gy  becomes 
0  <  gy  <8 \/jt2  ~  0.8106,  and  hence 

°dyn  <  l  +  — (for  appropriately  chosen  a ). 
^stat  ~  1  “  Pi 


There  is  a  large  difference  between  the  dynamic  stress  and  the  static  stress  as  the 
forcing  frequency  approaches  the  first  natural  frequency  (the  resonance  case).  Large 
differences  are  still  noticed  away  from  resonance.  For  example  at  /J,  =  0.5  the  dynamic 
stress  can  still  be  as  large  as  three  times  the  equivalent  static  stress  for  an 
inappropriately  chosen  strain  gauge  location. 

We  can  draw  several  conclusions  from  this  simple  example  comparing  the  stress  of  a 
dynamic  response  with  the  equivalent  static  case.  First,  if  the  vibration  originates  in 
another  part  of  the  truss,  then  the  only  way  it  can  enter  a  beam  under  consideration  is 
through  the  ends  (assuming  end  supports).  Assuming  vibrations  are  being  introduced 
from  both  ends,  then  Figure  C.2  suggests  the  best  place  for  the  strain  gauge  would  be 
in  the  middle.  However,  if  vibrations  are  introduced  only  from  one  end  of  the  beam, 
then  the  best  strain  gauge  location  (if  only  the  static  component  of  stress  is  required) 
would  be  at  the  opposite  end  of  the  beam.  The  above  analysis  also  suggests  that  the 
vibration  component  can  only  be  adequately  ignored  when  either  (or  both)  of  the  two 
following  conditions  is  satisfied: 

•  The  forcing  frequency  is  significantly  smaller  than  the  first  natural  frequency,  that 
is  fly  0.048  (or  /3y  ;S 0.075  if  a  is  appropriately  chosen).  (This  choice  of  yields 
a  maximum  error  of  10%  in  making  a  static  assumption  for  dynamic  loading). 

•  The  forcing  amplitude  is  smaller  than  the  static  loading  and  the  forcing  frequency 
is  smaller  than  the  first  natural  frequency. 
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D.  Least  Squares  Solutions 

In  this  appendix  we  review  some  least  squares  and  singular  value  decomposition 
theory. 


D.l  Singular  Value  Decomposition 

Following  the  outline  in  Golub  and  van  Loan  [33,  p.  71],  the  singular  value 
decomposition  (SVD)  of  a  real  m-by-n  matrix  A  is 

UtAV  =  D, 


where  the  matrix  D  is  diagonal  and  the  matrices  U  and  V  are  orthogonal,  and  have 
the  following  structure 


and 


D  -  diag(cr 


1  /•■■/Cpj 

P  = 

n 

3 

ll 

"o'1 

>» 

9T". 

min  \m,n}, 


The  function  "min {m,n}"  means  the  smaller  of  m  or  n.  The  singular  values  tr,  are 
ordered  so  that  ox  >o2  >...><7p  >0.  The  span (ulf...,ur)  and  span(ur+1,...,u„)  define 
the  range  and  nullspace  of  A  respectively,  where  r  =  rank(A)  (that  is, 
<7r+1  =  •••  =  a  =0).  The  2-norm  condition  number  of  the  matrix  A,  defined  in  terms  of  the 
singular  values,  is  given  by 


K 


that  is,  the  ratio  of  the  largest  to  smallest  singular  values.  The  2-norm  of  the  matrix  A  is 
simply  given  by  the  largest  singular  value,  that  is  ||A||2=<t1.  The  square  of  the 
Frobenius-norm  is  given  by  the  sum  of  the  singular  values  squared,  that  is 

l|A||p  =<71  +°2  +"'  +  <Jp  • 

Figure  D.l  shows  a  geometric  interpretation  of  singular  values  and  vectors  for  the  two- 
dimensional  case.  The  left  figure  shows  the  unit  circle,  which  represents  all  vectors  of 
unit  length.  While  the  right  figure  shows  the  mapping  of  the  unit  circle  under  the 
matrix  A  (a  linear  transformation).  The  singular  values  control  the  length  of  the  ellipses 
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axes.  In  general  a  singular  value  decomposition  may  be  thought  of  as  a  linear  mapping 
from  a  hyper-sphere  to  a  hyper-ellipsoid,  the  singular  values  controlling  the  length  and 
the  left  singular  vectors  controlling  the  direction  of  the  hyper-ellipsoid  axes.  Finally,  if 
the  left  singular  vectors  control  the  hyper-sphere's  elongation  directions,  then  the  right 
singular  vectors  control  the  remaining  transformations  of  the  sphere  (for  example 
rotations  and  reflections).  The  dashed  lines  show  the  effects  of  the  right  singular 
vectors. 


Figure  D.l.  Illustration  of  singular  values  and  vectors  for  the  two-dimensional 

case. 


If  r  is  the  rank  of  A  then  the  least  squares  (LS)  solution  is  given  by 


ujb 


x^=^-jrLv:=A+b 

1=1 


where  A+  is  the  pseudo-inverse  (defined  below).  The  solution  xLS  minimises  the 
2-norm  residual  ||Ax-b||2,  and  the  minimum  is  given  by 


m 

IIa*ls  -H2  =  ^{uJbf  ■ 

!=r+ 1 


The  pseudo-inverse  is  defined  as 


where 


A+=VD+Ut  , 


r 

D+  =  diag 

V 


.,0 


e  91” 
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Golub  and  van  Loan  [33]  use  the  pseudo-inverse  to  show  that  small  changes  in  A  or  b 
can  induce  arbitrarily  large  changes  in  the  least  squares  solution.  One  approach  to 
determine  the  numerical  rank  of  A  is  to  assign  a  tolerance  8  such  that  the  singular 
values  are  split  into  two  sets 

o1>...>or>8>or+1  >...>an, 

and  if  the  matrix  A  has  sd  significant  digits  one  suggested  choice  for  the  tolerance  is 
<5  =  l(Tsrf||A|L. 


D.2  Column  Weighting 

Since  the  loads  that  feed  into  the  system  are  of  different  natures  (for  example  bending 
moments,  shear  forces,  and  accelerations),  the  coefficients  of  the  resulting  stress 
equation  are  most  likely  to  be  different  orders  of  magnitude.  This  variation  in 
magnitude  leads  to  near  rank  deficient  load  matrices. 

Golub  and  van  Loan  [33,  p.  251]  state  that  the  least  squares  (LS)  solution  of  the  problem 
min||Ax  -  b||2 ,  where  Ae  9T'X"  and  b  e  9L" 

can  be  obtained  by  finding  the  minimum  2-norm  solution  yLS  to 

min||(AG)y  -  b\\2 

then  setting  xG  =  GyLS ,  where  G  e  9TX"  is  non-singular.  One  choice  of  G  which 
normalises  the  columns  is 


G  =  G0=diag(l/||A,1||2 . l/||Aj2), 

where  the  notation  Aij  denotes  the  zth  row  and  jth  column  of  the  matrix  A ,  and  thus 
G0  transforms  each  column  in  A  to  a  unit  vector  (under  a  2-norm).  However,  since 
column  weighting  affects  singular  values,  a  scheme  for  determining  numerical  rank 
may  not  return  the  same  estimates  when  applied  to  A  and  AG . 
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E.  Maximum  Error  of  an  Approximate  Matrix 


In  this  appendix  we  show  that  the  normalised  error  of  any  matrix  approximation  has 
two  minima  and  two  maxima.  Furthermore,  the  two  maxima  are  of  equal  magnitude 
but  point  in  exactly  opposite  directions,  with  a  similar  result  holding  for  the  two 
minima. 


Let  the  matrix  A  be  an  approximation  to  a  true  matrix  A .  Define  the  normalised  error 
squared  between  these  two  matrices,  as  a  function  of  some  vector  x ,  as 


e 


2 


(E.l) 


where  the  vectors  y  and  y  are  the  linear  transformations  of  the  vector  x  under  the 
matrices  A  and  A  respectively  (that  is  y  =  Ax  and  y  =  Ax).  Defining  the  error  matrix 
as 


E  =  A-A, 


then  Equation  (E.l)  may  be  re-written  as 


e 


2 


(E.2) 


Remembering  the  2-norm  relation  || zf2  =  zT  -z  we  can  easily  derive  the  quadratic  form 
(of  AtA) 

| \\Ax\f2  =  (Ax)r  -  (Ax) 

=  xtAtAx. 


Thus  we  see  that  the  normalised  error  (E.2)  will  have  exactly  two  maxima  and  two 
minima.  We  now  prove  this  result  for  the  special  case  of  a  two-by-two  matrix. 


Let  the  matrix  A  and  the  error  matrix  E  be  two-by-two  matrices  and  x  a  two- 
dimensional  vector,  each  containing  elements  a- ,  e(- ,  and  x,  respectively,  so  that 


flu 

U12 

,  E  = 

eu 

e12 

,  and 

x  = 

xa 

_a2\ 

a22_ 

e21 

e22_ 

_*2_ 

Then  the  quadratic  form  becomes 
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\\Axf2  =  xtAtAx 

=  (a21  +  ^21)^1  +2  (ana12  +  a21a71)x1x2  +  (a22  +a^)x \. 

Let  the  vector  x  be  a  unit  vector,  then  we  may  write  x  in  terms  of  a  real  parameter  6 
so  that 

x  =  [cos0  sin0]r. 


Note  we  have  not  lost  any  generality  in  setting  x  to  be  a  unit  vector  since  we  may 
divide  both  numerator  and  denominator  of  Equation  (E.2)  by  any  scalar.  Dividing  the 
numerator  and  denominator  of  the  normalised  error  (E.2)  by  the  scalar  ||x||2  forces  x  to 
be  a  unit  vector.  Using  the  above  unit  vector  definition  of  x  we  can  write  the 
normalised  error  in  terms  of  the  parameter  6 ,  which  gives 

g2  _  Of0  +  cq  cos 26  +  a2  sin  26  ^ 

P0  +  j8a  cos2 6  +  [52  sin2 6  ' 


where  the  constants  a,  are  a  function  of  the  error  matrix  E  only  and  are  given  by 

(Xq  —  1  "t  ^12  ^21  ^ 22  t 

2  2  2  2 
OC-y  —  ^21  ^22  7 


and 


a2  ~  2(ene12  +  e21e  22)  • 


Similar  expressions  arise  for  the  /?,  constants,  except  the  etj  are  replace  by  al} . 


We  want  to  determine  all  stationary  points  of  the  normalised  error  given  by 
Equation  (E.3),  so  as  to  locate  the  maximum  error.  We  assume  that  the  matrix  A  has  full 
rank  so  that  the  denominator  of  Equation  (E.3)  is  non-zero.  Differentiating  the 
normalised  error  with  respect  to  6  and  equating  the  result  to  zero,  after  a  small 
amount  of  simplification  we  are  left  with  the  condition 

y0  +  y1  cos  26  +  y2  sin  26  =  0 ,  (E.4) 


where 

7o  =aiPi  ~aiP2  /  7i  =«2Aj  and  72  =a0Ai  -«iA) 

are  constant.  We  will  show  later  that  we  need  only  consider  the  range  6  e  [0,;r) ,  so  that 
the  above  condition  has  only  two  roots  (one  points  in  the  direction  of  the  maximum 
and  the  other  in  the  direction  of  the  minimum).  The  condition  for  stationary  points  of 
the  normalised  error  given  by  Equation  (E.4)  can  be  shown  to  be  a  quartic  in  cos  2 6 
with  two  roots  (that  are  repeated)  given  by 


60 


DSTO-RR-0171 


COS  29  = 


-Y0Y1  ±Y2y[W+y[-ri 


2  2 

Yi+Yi 


For  the  denominator  in  the  above  expression  to  be  zero,  we  require  that  y\  +  y\  =  0 , 
which  implies  that  a0/f30  =  a1/p1  =  a2/p2 .  Thus  the  only  way  this  denominator  can  be 
zero  is  if  the  matrix  E  is  a  scalar  multiple  of  the  matrix  A,  in  which  case  the  normalised 
error  is  identically  equal  to  that  scalar  multiple  for  all  values  of  8  . 

Using  Equation  (E.2)  we  see  that  if  a  vector  x  maximises  the  error  giving  the  value  emax , 
then  the  diametrically  opposite  vector  -x  must  also  maximise  the  error  with  the  same 
maximum  value  emax .  This  justifies  the  earlier  restriction  of  the  parameter  6  to  the 
range  6e  [0,7r).  Using  the  fact  that  only  one  maximum  occurs  for  Be  [0,7r)  together 
with  the  fact  that  the  second  maximum  is  diametrically  opposite  to  this  first  maximum 
proves  that  the  normalised  error  surface  between  a  two-dimensional  matrix  A  and  its 
approximation  contains  two  maxima.  A  similar  result  holds  true  for  the  two  minima  of 
the  normalised  error  surface. 

We  can  generalise  the  above  result  for  the  n  -dimensional  case  by  holding  n  —  2 
variables  fixed  and  using  the  above  result  to  show  that  for  the  two  free  variables  the 
normalised  error  function  has  two  minima  and  two  maxima.  Any  search  algorithm 
then  proves  the  result  to  hold  for  all  higher  dimensions. 

Figure  2.1  shows  a  three-dimensional  example  for  three  orthogonal  search  paths.  We 
begin  by  choosing  an  arbitrary  great  circle  on  the  unit  sphere,  then  determining  the 
maximum  error  on  this  great  circle.  (This  initial  great  circle  and  corresponding  error 
function,  which  is  drawn  as  a  projection  orthogonal  to  the  sphere,  are  shown  in 
medium-grey.)  Note  that  the  maximum  was  simple  to  find  since,  as  we  proved,  the 
error  function  on  any  great  circle  has  only  two  maxima  (of  equal  magnitude  but 
diametrically  opposite  on  our  sphere).  Then  we  draw  a  second  great  circle  orthogonal 
to  the  first  great  circle  and  passing  through  this  error  maximum.  (This  second  great 
circle  and  corresponding  error  function  are  shown  in  light-grey.)  We  repeat  the 
previous  procedure,  finding  the  error  maximum  on  the  second  great  circle.  Then  draw 
a  third  great  circle  (shown  in  black)  orthogonal  to  the  second  great  circle  and  passing 
through  this  second  error  maximum.  Continuing  in  this  same  manner  we  eventually 
arrive  at  the  maximum.  A  simple  inductive  argument  carries  the  proof  to  all  higher 
dimensions. 

An  alternative  proof  is  to  consider  only  half  the  sphere.  Then  if  the  error  line 
corresponding  to  all  great  circles  on  this  hemi-sphere  have  only  one  maximum  (which 
we  have  already  proved),  the  error  surface  on  this  hemi-sphere  must  have  only  one 
maximum.  Again,  a  simple  inductive  proof  carries  results  onto  all  higher  dimensions. 
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Figure  El  Three-dimensional  visualisation  of  three  search  paths  (great  circles)  and 
their  corresponding  error  functions  (drawn  as  projections  emanating  orthogonally 
from  the  sphere).  The  first,  second,  and  third  search  paths  are  shown  in  medium- 
grey,  light-grey,  and  black  respectively. 


Consideration  of  the  geometric  interpretation  of  singular  values  (see  Appendix  D.l) 
shows  that  the  two  maxima  for  the  normalised  error  must  be  parallel  for  the  n- 
dimensional  case.  If  the  vector  location  of  one  of  the  maxima  is  known,  then  the  second 
maximum  is  simply  given  by  the  same  vector  in  the  reverse  direction.  Figure  E.2 
illustrates  this  point  for  the  two-dimensional  case. 


Figure  E.2.  Linear  transformation  of  unit  vectors  by  the  error  matrix  E  and  the 
true  matrix  A.  The  maximum  normalised  errors  lie  on  the  same  line  (radiating 
from  the  origin)  but  in  the  opposite  direction. 

Finally,  it  appears  possible  that  using  singular  value  decomposition  information  alone 
the  maxima  and  minima  of  the  error  surface  could  be  directly  found  without  resorting 
to  an  iterative  minimisation  or  maximisation  algorithm.  However,  this  approach  was 
not  investigated. 


62 


DSTO-RR-0171 


E.l  Regions  of  Maximum  Error 


In  the  preceding  section  we  showed  that  the  maximum  error  occurred  on  the  nth 
dimensional  hyper-sphere  given  by  the  unit  force  vector,  and  that  the  error  surface  was 
quartic  (with  repeated  roots).  Figure  E.3  shows  four  random  examples  of  the  quartic 
error  surface  for  the  two-dimensional  case.  The  inner  curve  is  the  unit  circle,  while  the 
outer  curve  is  the  quartic  error  surface. 


1 


2 


3 


4 


Figure  E.3  The  error  surface  corresponding  to  four  random  E  and  A  matrices. 
Ill-conditioned  A  matrices  usually  lead  to  localised  errors. 


The  condition  numbers  for  the  error  and  exact  matrices,  k(E )  and  k(A),  along  with  the 
error  bound,  <Tmax(E)/<7min(A),  for  the  four  examples  in  Figure  E.3  are  shown  in 
Table  E.l.  We  see  from  this  small  sample  that  it  is  the  condition  number  of  the  exact 
matrix  A  that  determines  whether  or  not  the  maximum  error  is  localised.  Thus  if  the 
exact  matrix  is  ill-conditioned  and  the  error  matrices  E  are  uniformly  distributed  and 
not  too  ill-conditioned,  then  on  the  average  we  can  expect  localisation  of  the  two 
maximum  error  regions. 


Example 

1 

2 

3 

4 

k(E) 

2.29 

52.1 

1.55 

1.34 

k(A) 

1.70 

6.41 

54.1 

2.97 

<7max(E)/tfmin04) 

0.210 

0.968 

5.06 

0.687 

Table  E.l  Condition  numbers  and  error  for  the  examples  in  Figure  E.3. 


Let  us  consider  several  different  scenarios  for  the  two-dimensional  case  to  further 
illuminate  the  above  statements.  Assume  that,  for  the  matrices  E  and  A,  all  these 
scenarios  have  the  same  set  of  singular  values  (SVs)  but  different  orientations  (refer  to 
Figure  D.l).  That  is,  in  the  terminology  of  Appendix  D.l,  the  SV  and  right  singular 
vectors  remain  the  same,  and  only  the  left  singular  vectors  change  between  scenarios. 
Furthermore,  so  that  the  matrix  A  dominates  the  response  of  the  error  surface,  assumed 
the  ratio  of  condition  numbers  (of  A  to  E)  is  large,  that  is  k(A)/k(E)»1.  (This 
assumption  corresponds  to  A  ill-conditioned  and  E  far  less  ill-conditioned.)  Then  the 
two  error  regions  surrounding  the  two  maxima  will  be  most  localised  when  the 
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orientations  of  ex,  (E)  (the  largest  SV  of  E)  and  <J2(A)  (smallest  SV  of  A)  are  collinear. 
At  the  other  extreme,  the  two  error  regions  surrounding  the  two  maxima  will  be  most 
smeared  when  the  orientations  of  <72(E)  (the  smallest  SV  of  E)  and  ct2(A)  (smallest  SV 
of  A)  are  collinear.  However,  it  is  the  matrix  A  that  dominates  the  control  of  the  error 
surface.  Hence  for  most  cases  it  is  the  amount  of  ill-conditioning  within  matrix  A  that 
will  affect  whether  or  not  the  region  of  maximum  error  is  localised.  Under  the  above 
assumptions,  we  also  expect  that  these  localised  error  regions  will  become  more 
pronounced  with  increasing  matrix  dimension. 
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